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1.  SCIENTIFIC  BACKGROUND 


1.1  Motivation 

Unlike  the  H 3  system  in  its  ground  electronic  state,  which  does  not  support  any  bound 
states  of  nuclear  motion,  this  system  is  known  to  have  a  long-lived1-4  excited  Rydberg  state 
2pz  2A-2.  The  high  resolution  spectroscopic  measurements  of  Herzberg  and  co-workers5-9 
and  of  Helm  and  co-workers10-13  confirmed  it  beyond  doubt.  A  very  intense  beam  of  Hz  in 
this  metastable  state  has  been  generated  by  Garvey  and  Kuppermann14.  Their  estimation 
of  the  lifetime  of  this  metastable  state  is  more  than  40  p  seconds.  This  metastable  Hz 
species  can  liberate  180  Kcal/mole  (7.8  eV/molecule)  when  decomposing  into  H^[X  l£+) 
molecules  according  to  HzfiPz  2A%)  — *  2/2H2{X  L££)-  Its  specific  impulse  (ISP)  is 
estimated  to  be  about  2050s,  while  current  rocket  fuels  have  ISP  of  the  order  of  only  400s. 
Such  a  high  ISP  makes  Hz  a  very  interesting  rocket  propellant  candidate. 

What  the  lifetime  of  this  metastable  Hz  actually  is,  via  what  kind  of  mechanism  it 
decays,  and  how  collisions  affect  this  lifetime  are  important  questions  whose  answers  are 
needed  in  assessing  its  potential  as  a  possible  rocket  propellant.  One  of  the  decay  channels 
is  the  ro-vibrational  predissociation  of  the  2pz  2A'£  metastable  electronic  state  into  the 
2p  2E'  repulsive  ground  electronic  state.  It  is  difficult  to  measure  directly  the  lifetime  of 
this  channel,  because  it  is  not  accompanied  by  any  emission  of  radiation.  The  total  lifetime 
estimation  from  linewidth  measurements5-9  suffers  from  the  Doppler  broadening  effects 
inside  the  plasma  sources  which  are  used  to  generate  Hz  metastable  molecules,  especially 
if  the  lifetime  of  interest  is  long,  as  is  the  case  for  the  Hz  in  the  2 pz  2A'£  state.  For  this 
reason,  a  theoretical  investigation  of  its  lifetime  (  for  both  radiative  and  predissociative 
processes  )  is  highly  desirable.  This  theoretical  investigation  should  also  provide  guidance 
for  the  experimental  studies  of  the  properties  of  this  metastable  species. 

1.2  Decay  Processes  and  Lifetime 

In  the  theory  of  predissociation  processes,  the  lifetime  r  of  a  bound  state  predissoci¬ 
ating  to  an  unbound  state  can  be  expressed  as15, 10 

C  =  A|  (1) 

Vg‘=<9m\H'\*ls>  (2)' 

In  these  expressions,  is  the  wavefunction  of  the  bound  state  with  quantum  number  m 
and  energy  f?m)  is  that  of  the  unbound  state  with  total  energy  E ,  and  H'  is  the  part 
of  total  Hamiltonian  which  provides  the  coupling  between  the  bound  state  and  unbound 
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state.  Vp  is  called  the  coupling  matrix  element  between  these  two  states.  The  integration 
in  Eq.  2  is  over  the  coordinates  of  all  electrons  and  nuclei.  The  choices  of  normalization 
for  the  above  two  wavefunctions  are: 


<  |  >=  1 


(3) 


<  |  *E'  >=6{E-E ') 


(4) 


For  electric  dipole  transitions,  the  Einstein  coefficient  between  an  initial  state  \Pt-  and 
a  final  state  /  is17 

Aiif  =  (4qw3/3c2)|<  'Hi  |  T  |  Vj  >|2  (5) 

and  the  spontaneous  radiation  lifetime  of  this  process  is 


r  = 


(6) 


where  a  is  the  fine  structure  constant,  w  the  frequency  of  the  emitted  photon,  and  c  the 
speed  of  light.  Here  again  the  initial  and  final  state  wavefucntions  are  normalized  as  in 
Eqs.  3  and  4.  <  |  T  |  Vf  >  is  the  electric  transition  dipole  moment  between  the 

initial  and  the  final  states,  and  again  the  integration  in  Eq.  5  is  over  the  coordinates  of  all 
electrons  and  nuclei.  For  molecular  systems,  the  total  wavefunctions  in  Eqs.  2  through  4 
can  be  written,  using  the  Bom-Oppenheimer  approximation18,  as 


'Ftata/  —  ^ elcctr^ nucl  (7) 

Here  '& rjectr  is  the  electronic  wavefunction  with  fixed  nuclear  coordinates  as  parameters 
and  ynuci  describes  the  nuclear  motion  in  an  effective  potential  generated  by  the  electron 
motion  (*.e.,  the  electronic  potential  energy  surface). 

In  order  to  calculate  both  the  radiative  and  predissociative  lifetimes,  Eqs.  1  through 
7  show  that  the  following  quantities  are  needed: 

1.  The  wavefunctions  which  are  solutions  of  the  electronic  motion  Schrodinger  equa¬ 
tion  for  the  2 pz  2 A %  state  and  all  lower  energy  electronic  states  (2p  2E'(  1), 
2p  2E'{ 2)  and  2s  2A\),  to  which  the  2 pz  2A%  state  may  decay  either  radiatively  or 
predissociatively.  The  resulting  electronic  potential  energy  surfaces  are  needed 
for  solving  for  the  wavefunctions  of  the  nuclear  motion. 
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2.  The  coupling  matrix  elements  between  those  electronic  states. 


3.  The  electric  transition  dipole  moments  among  those  electronic  states. 

4.  The  ro-vibrational  states  of  the  nuclear  motion  in  the  field  of  the  potential  energy 
surfaces  which  support  bound  states.  The  energy  eigenvalues  of  those  states  will 
determine  the  spectroscopic  transition  energies,  and  the  wavefunctions  are  needed 
in  the  lifetime  calculations. 

5.  The  scattering  wavefunction  on  the  ground  potential  energy  surface. 


The  following  sections  describe  in  detail  the  progress  achieved  so  far  in  obtaining 
these  different  quantities  necessary  for  the  calculation  of  the  radiative  and  predissociative 
lifetimes. 


2.  Calculation  of  the  Electronic  Wavefunctions  and  Energies  of  Hz 
2.1  General  Consideration 

The  first  step  toward  the  investigation  of  the  lifetime  of  the  2 pz  metastable  state 
of  Hz  is  to  calculate  the  electronic  potential  energy  surfaces  of  the  ground  and  the  first 
three  excited  electronic  states,  as  well  as  the  electronically  non-adiabatic  coupling  matrix 
elements  and  the  electric  dipole  transition  moments  between  those  electronic  states.  We 
originally  hoped  that  these  quantities  would  be  calculated  by  other  groups.  Since  this  has 
not  turned  out  to  be  the  case,  we  initiated  a  program  to  calculate  them  ourselves. 

Let  us  first  introduce  the  notation  for  those  states.  When  the  three  protons  form  an 
equilateral  triangle,  the  electronic  ground  state  is  degenerate  with  the  first  excited  state, 
and  the  two  generate  the  2p  2E'  representation  of  the  Dzh  point  group19.  The  second  and 
third  excited  electronic  states  are  classified  as  2s  2A\  and  2p2  2A'{.  The  x  and  y  axes  are 
in  the  molecular  plane,  while  the  z  axis  is  perpendicular  to  it.  (2s,  2 px,y,z)  are  the  state 
assignments  for  an  Li  atom  in  the  united-atom  limit  (UA).  We  still  use  this  notation  to 
identify  those  states  even  when  the  nuclear  geometry  is  no  longer  an  equilateral  triangle. 

The  potential  energy  surface  is  well  known  when  the  three  electrons  are  in  the  2 p  2E'{  1) 
ground  state19-22.  Using  a  functional  extrapolation  and  the  double  many  body  expansion 
method23-26,  Varanda  et  al.  also  obtained  the  potential  energy  surface  of  the  2p  2E'{2) 
state  which,  for  the  equilateral  triangle  nuclear  geometry,  is  degenerate  with  the  ground 
state20,22.  Along  with  the  potential  energy  surfaces,  the  major  coupling  matrix  elements 
between  the  ground  state  and  the  first  excited  state  were  also  obtained  using  the  same 
method22. 

After  the  first  Hz  spectroscopic  measurements3,  several  theoretical  electronic  calcu¬ 
lations  for  this  system  with  equilateral  triangle  nuclear  geometry  have  been  made27-30, 
Emphasis  was  placed  on  the  calculation  of  a  large  number  of  electronic  states,  in  order 
to  reproduce  the  observed  pattern  of  the  experimental  energy  levels  in  the  Hz  Rydberg 
spectra.  Figure  1  is  the  state  correlation  diagram  between  an  equilateral  triangular  Hz  and 
a  diatom  #2  plus  a  distant  atom  H.  In  these  calculations,  a  frozen  core  approximation 
was  used  for  H^  ,  plus  single  Rydberg  excitations  for  the  third  electron.  These  limited 
calculations  agree  with  each  other  and  are  qualitatively  successful  in  explaining  the  exper¬ 
imental  spectra.  The  frozen  core  treatment  gives  however  poor  results  for  the  2 p  2E'{  1,2) 
states,  and  also  over-estimates  the  energy  gap  between  the  2s  2A\  and  2 p2  2A'2  states  by 
a  factor  of  1.6  to  2.0  32  at  the  inter-nuclear  distance  R  =  1.64  bohr. 
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H.fX^  +  hT  +  e- 
H2  (x  'Zg) +  H(n-3) 

H2(X  'L*)  +  H  +  (n=2) 

3  H  (n=1) 

H2(X  ’l)  +  H(n=1) 


Figure  1:  Energy  Level  and  Correlation  Diagram  for  H3.  The  spacing  of  the  H3 
energy  levels  was  obtained  theoretically  for  an  equilateral  triangle  configuration27  and 
referred  to  the  energy  of  dissociated  products  by  the  results  of  a  separate  calculation31. 
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Recently  another  calculation  for  several  of  the  low  electronic  states  of  II3  in  some 
special  geometries  has  been  published32.  The  MRD-CI  method  used  in  that  work  is  very 
general  and  powerful33,34.  It  could  be  used  with  arbitrary  nuclear  geometries  and  to 
generate  the  full  potential  energy  surfaces.  It  also  offers  a  much  better  description  of  the 
2 p  2£'(l,  2)  states,  and  reduces  the  energy  gap  between  the  2s  2A\  and  2 pz  2A!\  states 
to  just  10%  greater  than  the  experimental  measurement32.  The  electric  dipole  transition 
moment  between  these  two  states  for  an  equilateral  triangle  nuclear  geometry  with  the 
intemuclear  distance  of  R=1.64  bohr  agrees  with  the  previous  calculations.  Because  of 
its  capability  and  performance,  we  chose  this  method  for  our  calculations  and  initiated  a 
collaboration  with  Professor  J.S.  Wright32  of  Carleton  University,  Ottawa,  Canada. 

2.2  MRD-CI  method 

The  details  behind  the  MRD-CI  (Multi-Reference,  single  and  Double  excitation  Con¬ 
figuration  Interaction)  are  beyond  the  scope  of  this  report  and  can  be  found  elsewhere33,34. 
Only  the  major  steps  of  the  calculation  are  briefly  described  below: 

1.  First,  a  set  of  atomic  orbitals  (AO)  are  chosen.  Gaussian-Type  atomic  orbitals 
(GTO)  are  used35. 

2.  Second,  the  nuclear  geometry  is  taken  into  account.  A  set  of  symmetry-adapted 
functions  (SAF)  is  constructed  by  linear  combinations  of  the  set  of  atomic  orbitals. 
By  taking  advantage  of  the  nuclear  geometry  symmetry,  the  MRD-CI  calculation 
for  a  given  molecular  system  can  be  speeded  up  significantly. 

3.  In  the  third  step,  the  Self  Consistant  Field  (SCF)  calculation  with  this  SAF- 
AO  basis  set  is  conducted  in  an  iterative  manner.  The  molecular  orbitals  (MO) 
obtained  from  the  SCF  serve  as  the  starting  point  for  the  Cl  calculation. 

4.  In  the  forth  step,  before  the  Cl  calculation,  a  reference  configuration  set  and 
an  threshold  energy  are  chosen.  All  the  single  and  double  excitations  over  the 
reference  configurations  are  generated.  First,  a  small  scale  Cl  calculation  with 
iVref  reference  configurations  is  performed  and  an  estimation  of  the  eigen-energies 
of  these  wanted  states  is  obtained.  Then,  one  by  one,  each  generated  configuration 
is  tested  by  being  added  to  the  reference  configuration  set  and  another  small  scale 
Cl  calculation  of  NTe{  +  1  configurations  is  performed.  Only  those  configurations 
which  lower  the  energy  of  any  one  of  the  states  of  interest  by  an  amount  bigger 
than  the  threshold  energy  are  included  in  the  final  Cl  calculation.  This  is  a 
way  of  drastically  reducing  the  final  Cl  space  size  without  missing  important 
contributions  from  the  full  single  and  double  excitation  Cl  space.  In  order  to 
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ensure  the  convergence  of  the  calculation,  the  reference  configuration  set  is  chosen 
to  be  big  enough  so  that  the  final  wavefunction  of  each  wanted  state  has  at  lease 
90%  of  its  contribution  from  the  reference  configurations. 

5.  The  fifth  step  uses  a  large  fraction  of  the  total  CPU  time.  It  is  the  diagonalization 
of  the  electronic  Hamiltonian  including  all  the  configurations  selected.  After 
that,  another  calculation  with  a  new  threshold  twice  as  large  as  old  one  is  made, 
including  the  selection  of  the  configurations  to  be  added  to  the  reference  set  and 
the  final  diagonolization.  Assuming  the  correctness  of  the  extrapolation  back 
to  the  zero  threshold  energy,  the  final  Cl  energies  can  be  approximated  to  even 
better  values,  without  doing  the  calculation  at  such  a  small  threshold  energy  as 
to  make  the  calculation  too  costly.  This  is  the  most  important  aspect  of  this 
MRD-CI  method. 

6.  With  the  resulting  electronic  wavefunctions,  the  electric  transition  dipole  moment 
between  any  two  known  electronic  states  is  then  calculated. 

2.3  Numerical  Details  and  Results 

The  version  of  the  MRD-CI  code  for  Sun  workstations  was  obtained  from  Professor 
Wright32.  The  CRAY  version  was  obtained  directly  from  Professor  Buenker's  group33,34. 
These  codes  do  not  have  the  capability  of  calculating  the  electronically  non-adiabatic 
coupling  matrix  elements  between  two  electronic  states,  and  will  have  to  be  modified  in 
the  future  to  permit  such  calculations,  which  is  needed  in  the  predissociative  lifetime 
calculation. 

Initially  we  made  small  testing  runs  on  Sun  workstations  (Sun-386i,  Sun-3,  and  Sun- 
4)  and  a  micro  VAX.  The  production  runs  have  been  done  on  the  JPL  CRAY-XMP  and 
the  SDSC  CRAY-YMP.  The  typical  CPU  time  on  the  CRAY-YMP  is  about  400s  for  a 
complete  calculation  at  a  single  nuclear  configuration,  using  the  largest  one  of  the  basis 
sets  described  below. 

The  details  of  the  Gaussian-Type  atomic  orbitals  (GTO)  are  given  in  Ref.  31.  In 
summary,  they  are  defined  as: 

=  Nngti  r’*’-1exp(— &r2)yi,m(0,<£)  (8)' 

where  Nngti  is  a  normalization  factor  and  the  spherical  harmonic  function, 

(r,  are  the  spherical  coordinates  of  a  point  with  respect  to  the  origin  of  the  GTO. 
ng  can  have  the  values  1,  2,  3,  ...  ,  l  can  have  the  values  0,  1,  2,  ...  (which  corresponds  to 
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s,  p,  d,  ...  orbitals),  and  m  can  vary  from  —l  to  l  in  steps  of  unity,  i  is  an  index  for  the 
exponent  &.  Sometimes  rn'}~lYittn.{0,<f>)  is  replaced  by  linear  combinations  of  terms  of  the 
form  xnzyn* zn* ,  where  (x,  y,  z)  are  Cartesian  coordinates  and  nx,nv,  nz  are  non-negative 
integers  whose  sum  equals  ng  —  1.  . 

In  some  applications  contracted  basis  functions  are  used.  They  are  defined  as  linear 
combinations  of  similar  GTO  functions  having  the  same  quantum  numbers  n0 ,  l,  m  but 
different  exponents  &. 

Xcontr  {jlgi  ■  i  ^  ^  X  {.ngi  fi >  (9) 

t 

The  Ci  are  called  the  coefficients  of  the  contraction. 

We  have  used  two  GTO  basis  sets.  The  first  one  is  a  variation  of  the  basis  set  used 
in  Ref.  28.  It  consists  of  (l0s/8s,  4 p)  orbitals,  in  which  there  are  eight  s-type  GTOs 
(one  of  them  being  a  contraction  of  three  s-type  GTOs),  amd  four  p-type  GTOs.  The 
second  basis  set  consists  of  (12s/7s,  4p,  Id)  orbitals,  with  seven  s-type  GTOs  (one  of 
them  contracted  from  six  s-type  GTOs),  four  p-type  GTOs  and  one  d-type  GTOs.  The 
main  part  of  it,  10s/5s,  4p,  Id,  is  adopted  from  one  used  in  a  Hf  calculation36,  with 
the  intention  to  provide  simultaneously  an  accurate  description  of  the  system’s  valence 
space  and  the  Rydberg  states  arising  from  n=2  hydrogen  orbitals.  In  order  to  improve  the 
description  of  the  Rydberg  states  of  the  Hz  system,  we  added  two  diffuse  Gaussian  s-type 
functions  into  that  basis  set,  with  exponents  0.01149  and  0.0042  (see  Tables  1  and  2  for 
details). 

Both  GTO  basis  sets  have  been  tested  and  calibrated  in  calculations  for  an  isolated  H 
atom,  and  an  isolated  H 2  molecule  with  bond  length  R=1.4  bohr.  The  results  are  shown 
in  Table  3.  The  H  atom  calculation  shows  that  both  sets  give  a  fairly  good  description 
of  the  £T(ls),  H(2s),  and  H{2px<,J<z)  atomic  states.  The  H2  results  clearly  show  that 
the  (12s/7s,  4p,  Id)  basis  set  offers  a  much  better  description  of  the  H2(X  *E+)  and 
H2(b  3E+)  states  than  does  the  (10s/8s,  4p)  one.  Because  we  are  interested  in  correlating 
the  2 pz  2^2)  2s  2A[  and  2p  2i?'(l,2)  states  of  H3  with  these  two  diatomic  states  mentioned 
above  plus  a  free  H  atomic  state  (see  Figure  1),  it  is  important  to  have  a  good  description 
of  these  diatomic  states. 

In  the  Hz  system,  the  most  general  geometric  symmetry  is  the  reflection  symmetry 
Cs  with  respect  to  the  plane  of  these  three  nuclei.  In  order  to  take  advantage  of  this 
symmetry,  we  placed  the  nuclei  in  the  x  y  plane,  with  one  nucleus  at  the  origin  of  the 
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Table  1:  (10s/8s,  4p)  Gaussian-type  basis  set. 


orbital 

i 

6 

Is 

i 

68.16 

0.002558 

2 

10.2465 

0.01938 

3 

2.34648 

0.0928 

2s 

1 

0.67332 

1.00000 

3s 

1 

0.22466 

1.00000 

4a 

1 

0.0822 

1.00000 

5s 

1 

0.0475 

1.00000 

6s 

1 

0.01875 

1.00000 

7s 

1 

0.0133 

1.00000 

8s 

1 

0.00525 

1.00000 

Ip 

1 

0.7 

1.00000 

2p 

1 

0.20 

1.00000 

3p 

1 

0.06 

1.00000 

4p 

1 

0.024 

1.00000 

Table  2:  (12s/7a  4 p,  Id)  Gaussian-type  basis  set. 


orbital 


Is 

1 

2 

3 

4 

5 

6 

837.22 

123.524 

27.7042 

7.82599 

2.6504 

0.938258 

0.000112 

0.000895 

0.004737 

0.019518 

0.065862 

0.178008 

2s 

1 

0.372145 

1.00000 

3s 

1 

0.155838 

1.00000 

4s 

1 

0.066180 

1.00000 

5s 

1  ‘ 

0.027580 

1.00000 

6s 

1 

0.011490 

1.00000 

7s 

1 

0.004200 

1.00000 

lp 

1 

1.6 

100000 

2p 

1 

0.40 

1.00000 

3p 

1 

0.09 

1.00000 

<p 

1 

0.025 

1.00000 

Id 

_ 1 

1 

r 

1.0 

1.00000 
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Table  3:  Calibration  of  the  basis  sets  for  isolated  //  and  H 2“. 


Calculated  H  energy  levels 


Gaussian 
basis  set 

H(  Is) 

H(  2s) 

H(  2Px) 

tf(2p„) 

H{  2Pi) 

10s/8s,  4p 

-0.499942 

-0.124989 

-0.124812 

-0.124812 

-0.124812 

12s/7s,  4p,  Id 

-0.499998 

-0.124992 

-0.124723 

-0.124723 

-0.124723 

Exact  value 

-0.500000 

-0.125000 

-0.125000 

-0.125000 

-0.125000 

Calculated  H 2  energy  levels  (R=1.4) 


Gaussian 

basis  set 

h2(X  l2+) 

H2{b  3E:) 

10s/8s,  4p 

-1.170045 

-0.782718 

12s/7s,  4p,  Id 

-1.173652 

-0.783904 

KW6 

-1.174474 

-0.784150 

a:  The  energies  are  in  hartree  and  the  bond  lenth  of  diatom  is  in  bohr. 
b:  Accurate  values  from  W.  Kolos  and  L.  Wolniewicz37. 
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(x,  y,  z)  coordinate  system,  another  along  the  positive  x  axis  and  the  third  in  the  x  y  half 
plane  with  y  positive.  The  ground  state  2 p  2E,{  1),  and  the  excited  states  2 p  ~E'( 2),  and 
2s  2A[  are  symmetric  under  that  symmetry  operation  (z  — ►  —  z),  while  the  2p-  -vlj  state 
is  anti-symmetric. 

We  place  sets  of  GTOs  (as  listed  in  Tables  1  and  2)  at  each  of  the  three  nuclei.  There 
are  20  basis  functions  for  the  (10s/8s,  4 p)  basis  set  and  25  for  the  (I2s/7s,  4p,  Id)  basis 
set,  where  there  is  one  function  for  each  s  type  GTO,  three  for  each  p  type  GTO  and  six 
for  each  d  type  GTO.  Because  we  have  three  nuclei,  the  total  numbers  of  basis  functions 
are  60  for  the  first  case  and  75  for  the  second. 

We  obtain  the  energies  and  wavefunctions  for  the  ensemble  of  states  symmetric  with 
respect  to  reflection  through  the  nuclear  plane  in  one  calculation  and  the  antisymmetric 
ensemble  in  another  one.  The  energy  threshold  is  2.0  x  10“6  hartree  for  the  symmetric 
ensemble,  for  both  sets  of  basis  functions.  For  the  antisymmetric  ensemble,  the  energy 
threshold  is  0.5  x  10-6  hartree  for  the  (12s/7s,  4 p,  Id)  basis  set  and  1.0  x  10-G  hartree  for 
the  (l0s/8<!,  4p)  one.  Since  the  choice  of  reference  configurations  depends  on  the  nuclear 
geometry,  we  had  to  adjust  it  through  a  trial  and  error  process  for  each  nuclear  geometry. 
Generally  speaking,  we  used  about  45  to  49  reference  configurations  in  the  calculation  of 
the  symmetric  ensemble,  and  about  20  in  the  calculation  of  the  antisymmetric  one. 

In  Tables  4  and  5,  we  list  the  energies  of  the  first  four  electronic  states  with  the 
equilateral  triangle  nuclear  geometry  for  each  basis  set.  These  results  are  plotted  in  Figures 
2  and  3.  In  Tables  6  and  7,  we  list  the  electric  transition  dipole  moments  for  this  geometry. 
Comparison  of  the  results  from  the  (I2s/7s,  4p,  Id)  basis  set  for  the  2p  2£7/ (1,2)  states 
with  those  of  the  DMBE  calculation22  shows  that  our  curve  is  nearly  parallel  to  the  DMBE 
one  with  an  upward  shift  of  about  0.003  hartree.  Comparison  between  the  results  of  the 
(10s/8s,  4p)  basis  set  and  those  in  Ref.  28  shows  good  agreement  (except  at  R  =  2.6 
bohr).  GTO  functions  located  at  the  the  center  of  mass  of  three  nuclei  have  been  used 
in  Ref.  28,  while  a  larger  set  of  basis  functions  located  at  each  nucleus  is  used  in  present 
calculation.  Comparison  of  Tables  4  and  5  shows  that  the  results  with  the  two  basis  sets 
listed  in  Tables  1  and  2  are  not  parallel  to  each  other.  This  indicates  that  it  is  necessary 
for  us  to  use  the  larger  basis  set  (I2s/7s,  4p.  Id)  in  our  fined  production  runs  in  order  to- 
get  a  potential  energy  surface  having  the  correct  shape. 

The  Cg  symmetry  ensures  that  the  electric  transition  dipole  moments  between  the 
antisymmetric  2 pz  2 A %  state  and  the  symmetric  2p  2E'{  1),  2p  2E'{ 2)  and  2s  2A\  states  have 
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Table  4 


Electronic  potential  energies  (in  hartree)  of  the  first  four  electronic  states  of  //3a  for  the  (l0s/8s,  4 p)  bas 


R6 

2p  2E'(  1) 

2p  2E( 2) 

2s  2A\ 

2 p,  7A'l 

1.0 

-1.273839 

-1.273825 

-1.267978 

-1.245835 

1.2 

-1.431537 

-1.431503 

-1.404972 

-1.388000 

1.4 

-1.510877 

-1.510886 

-1.462253 

-1.451243 

1.6 

-1.549290 

'  -1.549231 

-1.477925 

-1.470966 

1.64 

-1.553811 

-1.553849 

-1.478022 

-1.471764 

1.8 

-1.564752 

-1.564771 

-1.471242 

-1.467687 

2.0 

-1.568660 

-1.568710 

-1.452355 

-1.451797 

2.2 

-1.565200 

-1.565174 

-1.427605 

-1.420583 

2.4 

-1.558467 

-1.558474 

-1.400492 

-1.404303 

2.6 

-1.550885 

-1.550780 

-1.373003 

-1.378150 

2.8 

-1.542112 

-1.542476 

-1.345840 

-1.352258 

3.0 

-1.534483 

-1.534565 

-1.319509 

-1.327177 

a:  The  origin  of  energy  is  that  of  the  six  particles  (three  electrons  and  three  protons) 
at  infinite  separation.  The  energy  of  three  independent  f/ (Is)  is  -1.500000  hartree. 
These  three  nuclei  form  an  equilateral  triangle. 

b:  in  bohr. 
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Table  5 


Electronic  potential  energies  (in  hartree)  of  the  first  four  electronic  states 
of  Hza  for  the  (12s/7s,  4 p,  Id)  basis  set. 


R6 

2p  2E'(1) 

2 p  2E'( 2) 

2s  2A\ 

2 P,  2 A'' 

1.0 

-1.286448 

-1.286430 

-1.280663 

-1.258265 

1.2 

-1.441703 

-1.441650 

-1.415028 

-1.398848 

1.4 

-1.518046 

-1.518017 

-1.468988 

-1.458043 

1.6 

-1.554349 

-1.554268 

-1.482113 

-1.475586 

1.64 

-1.558556 

-1.558507 

-1.481895 

-1.475980 

1.8 

-1.569022 

-1.568989 

-1.474258 

-1.471001 

2.0 

-1.571945 

-1.571928 

-1.455205 

-1.454669 

2.2 

-1.568548 

-1.568561 

-1.430550 

-1.432079 

2.4 

-1.561349 

-1.561420 

-1.403023 

-1.406783 

2.6 

-1.552813 

-1.552907 

-1.375206 

-1.380527 

2.8 

-1.544312 

-1.544450 

-1.347990 

-1.354630 

3.0 

-1.536907 

-1.536859 

-1.322044 

-1.329407 

a:  The  origin  of  energy  is  that  of  the  six  particles  (three  electrons  and  three  protons) 
at  infinite  separation.  The  energy  of  three  independent  H (is)  is  -1.500000  hartree. 
These  three  nuclei  form  an  equilateral  triangle. 

b:  in  bohr. 
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ro 


R  (bohr) 


Figure  2:  Potential  Energy  Curves  for  Equilateral  H3.  R  is  the  length  of  the  side  of 

the  triangle.  The  (10s/8s,  4 p)  basis  set  was  used.  The  2 p  2i£'(l)  and  2 p  2E‘(2 )  curves 
are  not  completely  degenerate  with  each  other  because  of  inaccuracies  introduced  by  the 
limited  basis  set.  But  the  difference  between  them  can  not  be  discerned  on  the  scale  of, 
the  plot  (see  Table  4). 


15 


Figure  3:  Potential  Energy  Curves  for  Equilateral  H3.  R  is  the  length  of  the  side  of 
the  triangle.  The  (12s/7s,  4 p,  Id)  basis  set  was  used.  The  2 p  2.E'(l)  and  2 p  2E'( 2)  curves 
are  not  completely  degenerate  with  each  other  because  of  inaccuracies  introduced  by  the 
limited  basis  set  (see  Table  5). 


Table  6 


Transition  dipole  moment*1  (in  atomic  units)  among  the  first  four  electronic  states 
of  Hz  for  the  (10s/8a,  4p)  basis  set. 


R6 

T41(») 

T42(z) 

T43(z) 

1.0 

0.769(-3) 

0.4Sl(-3) 

2.61 

1.2 

0.230(-3) 

2.63 

1.4 

0.642(-4) 

0.705(-3) 

2.65 

1.6 

0.101(-2) 

-0.593(-4) 

2.66 

1.64 

-0.108(-5) 

0.148(-3) 

2.66 

1.8 

0.376(-3) 

-0.486(-3) 

2.66 

2.0 

-0.159(-3) 

0.204(-3) 

2.69 

2.2 

0.596(-3) 

-0.751(-4) 

2.70 

2.4 

0.120(-2) 

-0.116(-2) 

2.73 

2.6 

0.707 (-3) 

0.102(-2) 

2.76 

2.8 

0.531(-3) 

0.677(-3) 

2.78 

3.0 

-0.349(-3) 

0.27l(-3) 

-2.80 

Rh 

T31(x) 

T31(y) 

T32(x) 

T32(y) 

T21(x) 

T21(y) 

1.0 

1.70 

— 

1.48 

1.70 

-0.759(-2) 

-0.590(-l) 

1.2 

1.64 

I 

-1.24 

1.64 

-0.269(-l) 

0.935(-l) 

1.4 

-0.673(-l) 

1.85 

-0.651(-1) 

0.139 

0.974(-2) 

1.6 

1.44 

-0.803 

0.803 

1.44 

-0.938(-l) 

-0.152 

1.64 

1.11 

1.19 

-1.19 

1.11 

0.136(-l) 

0.186 

1.8 

-0.628 

1.35 

1.35 

0.627 

-0.142 

0.168 

2.0 

1.29 

-0.290 

.290 

1.29 

-0.108 

2.2 

1.12 

0.360 

-0.359 

1.12 

SH 

0.158 

2.4 

0.953 

0.443 

-0.440 

0.956 

-0.179 

0.212 

2.6 

0.373 

0.840 

0.841 

-0.371 

-0.180 

-0.198 

2.8 

0.393 

0.733 

0.730 

-0.389 

-0.140 

-0.209 

3.0 

-0.177 

-0.720 

-0.718 

0.174 

-0.194 

-0.101 

a:  Tij  is  the  transition  dipole  vector  for  equilateral  triangle  configurations.  The  indices  (1,  2,  3,  4)  refir 
to  the 

four  electronic  states  (2p  2£'(1),  2p  2E'{ 2),  2s  2A\,  2 pz  2 A 2)  respectively, 
b:  in  bohr. 
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Table  7 


Transition  dipole  moment"  (in  atomic  units)  among  the  first  four  electronic  states 
of  H3  for  the  (12a/7s,  4 p,  Id)  basis  set. 


R6 

T41(i) 

T42(z) 

T43(s) 

1.0 

-0.432(-3) 

-0.194  (-3) 

-2.61 

1.2 

0.541(-3) 

0.120(-3) 

-2.63 

1.4 

-0.483(-3) 

-0.489(-3) 

-2.65 

1.6 

0.809(-3) 

0.546(-3) 

-2.68 

1.64 

-0.830(-3) 

0.378(-3) 

2.69 

1.8 

0.647  (-3) 

-0.538(-3) 

-2.71 

2.0 

0.497(-3) 

0.732(-3) 

-2.74 

2.2 

-0.140(-3) 

0.100(-2) 

-2.75 

2.4 

0.904(-3) 

-0.556(-3) 

-2.78 

2.6 

0.152(-2) 

-0.153  (-3) 

-2.80 

2.8 

~0.206(-3) 

-2.81 

3.0 

0.710(-3) 

-0.758(-4) 

-2.82 

R6 

T31(x) 

T31(y) 

T32(x) 

T32(y) 

T21(x) 

T21(y) 

1.0 

-2.25 

-0.944(-l) 

-0.927(-l) 

2.26 

0.586(-l) 

-0.521(-1) 

1.2 

0.199 

2.04 

2.04 

-0.198 

-0.916(-1) 

-0. 182(- 1) 

1.4 

-0.161 

-1.83 

-1.84 

0.159 

-0.131 

-0.227(-l) 

1.6 

1.22 

1.10 

1.09 

1.22 

-0.202(-l) 

-0.172 

1.64 

-1.38 

0.801 

0.803 

1.39 

0.909(-l) 

0.157 

1.8 

0.362 

1.41 

1.41 

0.361 

0.189 

0.103 

2.0 

0.372 

1.25 

1.25 

-0.371 

-0.207 

-0.135 

2.2 

0.424 

1.07 

1.06 

-0.421 

-0.194 

-0.182 

2.4 

0.474 

0.918 

-0.917 

0.472 

0.161 

0.225 

2.6 

0.446 

0.812 

-0.814 

0.444 

0.148 

-0.226 

2.8 

0.443 

0.702 

0.700 

-0.439 

-0.111 

-0.226 

3.0 

-0.388 

-0.639 

-0.632 

0.385 

-0.104 

-0.196 

a:  Tij  is  the  transition  dipole  vector  for  equilateral  triangle  configurations.  The  indices  (1,  2,  3,  4)  refer 
to  the  four  electronic  states  (2p  2£'(1),  2p  2E'{ 2),  2s  2A\,  2 px  2A'£)  respectively. 

b:  in  borh. 
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only  z  components,  and  the  ones  between  these  symmetric  states  have  no  z  components. 
Since  the  wavefunctions  have  been  determined  by  the  variational  calculation  up  to  a  phase 
factor,  all  transition  dipoles  are  subject  to  a  possible  sign  change. 

Although  C a  is  the  only  symmetry  embedded  into  the  calculation,  when  three  nuclei 
form  an  equilateral  triangle,  the  Dzh  symmetry  group  associated  with  this  geometry  will 
manifest  itself  via  the  following  features: 

1:  The  2 p  2E'(  1)  and  2p  2E'{ 2)  states  are  nearly  degenerate. 

2:  The  electric  transition  dipoles  from  the  2 pz  2A %  state  to  the  2p  2E'(  1)  and  2 p  2E'( 2) 
states  are  close  to  zero. 

3:  Because  of  the  degeneracy  between  the  2p  2E’(  1)  and  2 p  2E'2  states  (under  the  sym¬ 
metry  of  an  equilateral  triangle),  they  can  always  be  written  as: 

|  2p  2£'(1))  =  cos ifi  |  <t> i)  -I- simp  |  fa)  (10) 

|  2 p  2E'{ 2))  =  -sinp  |  4>i)  +  cosp  \  <fo)  (11) 

|  4>\)i  |  <t> 2)  are  solutions  of  the  electronic  wave  equation  with  the  same  energy,  which 
form  another  E'  representation  of  the  Dzu.  group.  The  phase  p  is  not  determined  by 
the  variational  method  alone,  and  can  have  an  arbitrary  value.  For  two  calculations 
with  different  inter-nuclear  distances,  the  relative  phase  of  these  two  electronic  cal¬ 
culations  is  random,  which  in  turn  causes  the  x  and  y  components  of  the  transition 
dipole  moments  (T31,  T32,  and  T21)  to  vary  greatly  (see  Tables  6  and  7).  Even  so, 
the  symmetry  ensures  that: 

•  The  magnitudes  of  T31,  T32,  and  T21  do  not  depend  on  the  phase  p  and  thus 
change  smoothly  with  the  inter-nuclear  distance. 

•  i  T31  |  =  i  T32  I,  |  T31(x)  |  =  |  T32 (y)  j,  and  |  T31(y)  |  =  |  T32(x)  | 

All  of  these  features  have  been  confirmed  numerically  by  the  results  in  Tables  6  and  7 
and  by  Figures  4  and  5.  Since  the  molecular  properties  are  more  sensitive  to  the  quality  of 
the  wavefunction  than  the  energy  eigenvalues  are,  the  results  of  the  electric  transition  dipole 
moment  calculations  offer  another  strong  indication  that  the  wavefunctions  we  obtained 
are  of  good  quality. 

Our  results  for  the  electric  transition  dipole  between  the  2s  2A[  and  2 pz  2A2  states 
(T43(z))  show  a  very  slow,  smooth  and  monotonic  variation  with  the  size  of  the  equilateral 
triangle.  It  is  interesting  to  note  that  the  value  of  this  transition  dipole  moment  is  close 
to  that  between  the  2s  and  2 pz  states  of  an  isolated  H  atom  (  3.00  atomic  units).  The 
correlation  diagram  of  Figure  1  shows  that  those  two  states  dissociate  into  Hz(X  x£+)  + 
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Figure  4:  Length  of  the  transition  dipole  moment  T31  for  Equilateral  H3.  R  is  the 
length  of  the  side  of  the  triangle.  |T31|  is  the  length  of  the  T31  transition  dipole  moment 
(see  Tables  6-7  for  the  defination).  The  lower  curve  was  obtained  with  the  (I0s/8s,  4 p) 
basis  set  used,  and  the  upper  one  with  the  (12s/7s,  4p,  Id)  basis  set  used. 
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Figure  5:  Length  of  the  transition  dipole  moment  T21  for  Equilateral  H3.  R  is  the 
length  of  the  side  of  the  triangle.  |T21|  is  the  length  of  the  T21  transition  dipole  moment 
(see  Tables  6-7  for  the  defination).  The  lower  curve  was  obtained  with  the  (10s/8s,  4p) 
basis  set  used,  and  the  upper  one  with  the  (12s/7s,  4p,  Id)  basis  set  used. 
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//(2s)  and  H2[X  lE+)  +  H[2pz)  respectively,  which  should  indeed  furnish  transition 
dipole  moments  around  3.00  a.u.  As  a  result  we  can  expect  that  this  transition  dipole 
moment  will  be  more  or  less  the  same  for  all  the  nuclear  geometries  we  are  interested  in. 
This  turns  out  to  be  very  important  in  making  correct  state  assignments  in  our  calculation. 

At  the  equilateral  triangular  geometry  with  an  inter- nuclear  distance  R  =  1.64  bohr 
(  corresponding  approximately  to  the  equilibrium  geometry  of  the  metastable  i/3  ),  we 
calculated  the  energy  spacing  between  states  2s  2A[  and  2p  2 A %  to  be  1299  cm-1  with 
the  larger  basis  set  and  1374  cm-1  with  the  smaller  one,  while  the  best  value  previously 
calculated32  is  1422  cm-1  and  the  experimentally  estimated  value32  1256  cm-1.  Both  of 
our  results  for  the  electric  transition  dipole  moments  between  these  two  states  agree  with 
the  previous  calculation32  within  one  percent. 

In  Table  8,  we  present  the  results  of  of  the  energy  calculations  with  the  bond  angle 
fixed  at  60°  and  one  bond  length  fixed  at  10.0  bohr.  Since  the  H  atom  is  now  far  away 
from  the  Hi  diatomic  molecule,  we  expect  that  these  states  of  interest  will  correlate  with 
dissociated  states  according  to: 


H3{2p  2E'{  1))  —  H2{X  l£+'  //(Is) 

(12) 

#3 (2s  X)  — -  XL£)  +ZZ(2s) 

(13) 

— +Z/2(A’1E+)  +  Z/(  2Px) 

(14) 

—+»2{X  l;)  +  //(2Ps,) 

(15) 

zr3(2p  2a")  —  H2(X  XE+)  +  H(2Pz) 

(16) 

H3(2p  2E'( 2))  —  H2(b  3E^)  +  //(Is) 

(17) 

Since  the  energies  corresponding  to  the  right  side  of  Eqs  13,  14  and  15  are  the  same, 
we  write  Eqs  14  and  15  as  a  reminder  that  we  must  make  the  correct  state  assignments. 
The  energy  output  listed  in  Table  8  confirms  that  the  above  equations  describe  the  correct 
correlations  (also  see  Figure  6).  This  is  another  indication  that  our  GTO  basis  set  is  also 
appropriate  for  the  H2  molecule. 

In  the  equilateral  triangle  geometry,  the  energy  of  the  2 p  2E'{ 2)  state  is  below  that  of 
the  2s  2A\  state.  However,  for  configurations  in  which  H  and  H2  are  far  away  from  each 
other,  this  is  not  generally  true  (see  Figures  2,  3  and  6).  The  crossing  of  the  corresponding 
potential  energy  surfaces  increases  the  complexity  of  our  calculation,  because  we  have  to 
be  able  to  identify  the  states  correctly,  both  for  state  assignment  purposes  and  the  correct 
calculation  of  the  electric  transition  dipoles. 
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Table  8 


Potential  energies  when  H  and  Hn  are  far  away  from  each  othera. 


2p  2E'(l) 

2 p  2E'( 2) 

2  Pm 

2^' 

R“ 

this  work 

KWfc 

this  work 

KWC 

this  work 

KW'1 

1.2 

-1.663260 

-1.664934 

-1.218675 

-1.218964 

-1.288093 

-1.298934 

1.4 

-1.673015 

-1.674474 

-1.283682 

-1.284150 

-1.297616 

-1.299474 

1.6 

-1.667196 

-1.668580 

-1.331255 

-1.331724 

-1.292109 

-1.293580 

1.8 

-1.653721 

-1.655067 

-1.367693 

-1.368291 

-1.278657 

-1.280067 

2.0 

-1.636826 

-1.638132 

-1.396762 

-1.397064 

-1.261499 

-1.263131 

3.0 

-1.556159 

-1.557312 

-1.471696 

-1.472010 

-1.181171 

-1.182312 

a  :  One  bond  length  ia  fixed  at  10.0  Bohr.  The  bond  angle  is  60.0  degree. 

R  is  the  bond  length  of  /f2  (in  bohr). 

b  :  Potential  energy  of  H?(X  1S^')  from  W.  Kolos  and  L.  Wolniewicz37  with  the  energy  of  H(n  =  1) 
added. 

c  :  Potential  energy  of  ff^lb  3E„  )  from  W.  Kolos  and  L.  Wolniewic*37  with  the  energy  of  H(n  =  1)  added. 

d  :  Potential  energy  of  H?(X  iE+)  from  W.  Kolos  and  L.  WoLniewicz37  with  the  energy  of  H[n  =  2) 
added. 
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Figure  6:  Potential  Energy  Curves  for  Asymmetric  H3.  The  geometry  is  defined  in 
Table  8.  R  is  the  length  of  the  smallest  side  of  the  triangle.  The  largest  side  of  the  triangle 
is  10.0  bohr  and  the  bond  angle  between  them  is  60.0°.  The  (l2s/7s,  4p,  Id)  basis  set 
was  used. 
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In  order  to  map  out  the  interesting  part  of  the  potential  energy  surfaces,  wc  estimate 
that  about  500  nuclear  configurations  with  different  shapes  and  sizes  should  be  treated. 
We  have  at  this  moment  about  95%  of  them  done.  After  these  calculations  are  finished, 
they  will  be  fitted  by  a  procedure  developed  previously38  before  being  used  in  the  bound 
and  scattering  states  calculation  of  the  nuclear  motion. 

The  comparison  between  our  ab  initio  results  and  these  results  from  the  functional 
extrapolation  and  double  many  body  expansion  (DMBE)22  will  permit  us  to  check  if 
this  extrapolation  is  valid  for  nuclear  configurations  appreciably  removed  from  equilateral 
triangular  nuclear  geometries.  Recently,  the  H  +  H?  system  has  been  the  subject  of 
transition  state  spectroscopy  (TSS)  studies30-42.  There  are  some  theoretical  calculations 
on  the  spectroscopy  of  the  transition  between  the  ground  electronic  state  2 p  2E'{  1)  and 
the  excited  electronic  state  2 pz  but  they  are  limited  to  either  linear  configurations  or 

to  3 D  calculations  using  an  ad  hoc  potential  surface  for  the  excited  state  and  a  constant 
electric  transition  dipole  moments43.  Our  calculations  of  the  excited  surfaces  and  the 
nuclear  geometry-dependent  electric  transition  dipole  moments  between  them  will  be  very 
helpful  for  such  theoretical  studies  of  the  transition  state  spectroscopy  of  the  H  +  J?2 
system. 
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3.  Ro- vibrational  Bound  State  Calculation  for  the 
Upper  Manifold  of  the  Electronic  2 p  2E'  States  of  H 3 

The  next  step  towards  the  lifetime  calculation  is  to  obtain  the  ro-vibrational  nuclear 
wavefunctions  on  the  initial  and  final  electronic  potential  surfaces.  The  2 pz  2A,{  state 
of  H3  supports  bound  ro-vibrational  states  which  must  be  calculated.  There  are  several 
methods  for  performing  such  calculations44-48,  most  of  which  are  based  on  the  variational 
principle.  A  method  of  choice  should  have  the  following  properties: 

•  It  should  be  applicable  to  all  nuclear  geometries. 

•  It  should  be  able  to  treat  large  vibrational  amplitudes. 

•  It  should  be  able  to  take  advantage  of  special  nuclear  permutation  symmetries  of  the 

triatomic  system. 

•  It  should  be  robust,  and  easy  to  use. 

After  a  survey  of  the  available  methods,  we  chose  the  one  developed  by  Tennyson 
and  Sutcliffe45.  Since  the  only  previously  calculated  potential  energy  surface  of  H3  which 
supports  bound  ro-vibrational  states  (neglecting  its  coupling  to  the  ground  state)  is  the 
upper  manifoid  of  the  DMBE  2p  2E'  surface22,  we  used  it  as  a  testing  ground  for  the  bound 
ro-vibrational  state  calculation.  Because  the  conical  intersection  between  this  surface 
and  the  ground  state  one  introduces  the  Molecular  Aharonov-Bohm  (MBA)  effect49-52 
which  needs  special  treatment,  and  since  we  also  want  to  embed  the  correct  P3  nuclear 
permutation  symmetry  into  the  ro-vibrational  wavefunction,  we  also  developed  a  new 
hyperspherical  coordinate  propagation  method  which  is  very  general  and  powerful.  After 
obtaining  the  2s  2A[  and  2 pz  2A'2  potential  energy  surfaces,  we  will  use  those  two  methods 
to  calculate  the  corresponding  ro-vibrational  states. 

3.1  Variational  Method 

The  coordinates  used  in  the  method  of  Tennyson  and  Sutcliffe  are  the  scattering 
coordinates45.  Let  Ai,  A?,  A3  be  the  three  atoms  of  the  system,  and  (A ,»/,«)  be  any  cyclic 
permutation  of  (  1,  2,  3  ).  r>  is  the  internuclear  vector  for  the  diatom  AUAK  and  R*  the 
vector  of  A\  with  respect  to  the  center  of  mass  of  AUAK.  After  separating  the  motion 
of  the  center  of  mass  of  these  three  atoms,  the  system  can  be  described  by  six  variables 
(<*a,  0x,  lx,  R\,  rx,  0a).  The  three  angles  {ax,0x,lx)  are  Euler  angles,  which  describe 
the  orientation  of  the  tratomic  molecule  in  space.  r>  and  Ra  are  the  lengths  of  vectors  Ta 
and  Ra  respectively  and  Ox  is  the  angle  between  these  two  vectors. 

A  suitable  symmetrized  angular  basis  set  for  the  variational  calculation  is  chosen  to 
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be: 

|j,  k)  =(l+4,o)"'/2  2-  l/" 

(ei*Wi)i,iK,fc,1>)  +  (-l)P9y.-ib(«A)I>i,-k(a»,/3A,lA)}  (18) 

where  £fci0  equals  1  when  *  =  0  and  zero  otherwise.  DJM  *(«*>,/?>,  'lx,)  is  the  Wigner 
rotation  function53  and0y,*(0,x)  is  the  associated  Legendre  function54,  p  is  a  quantum 
number  that  can  assume  the  values  0  or  1.  J  is  the  total  angular  momentum  quantum 
number,  with  M  and  k  being  the  quantum  numbers  of  its  projections  along  the  space-fixed 
z  axis  and  the  body-fixed  z  axis  respectively,  j  is  the  angular  momentum  quantum  number 
of  the  diatomic  vector  rx .  The  Euler  angules  are  chosen  in  such  a  way  that  the  body-fixed 
z  axis  is  along  the  R>  direction  and  has  a  positive  projection  on  the  body-fixed  i  axis, 
(p,  J,  M)  are  constants  of  the  motion  for  the  system.  The  allowed  values  of  j  and  k  are: 

*=  (-7,-7+ 1,...,J- 1,J)  (19) 

i=(|*|,|*|  +  l,|*|  +  2,...)  (20) 


The  total  parity  of  the  spacial  wavefunction  under  inversion  through  the  system’s  center 
of  mass  is  (— 1) J+p. 

The  basis  functions  for  (Rx,  rx)  are  chosen  to  be  the  analytic  Morse  oscillator-like 
functions: 


where 


^,„(r>,  Rx)  =  j^H^H^Rx) 

(21) 

m,n  =  {0,1, 2, 3, 4,...} 

(22) 

Hn(r)  =  /?1/2lVn,aexp(-y/2)yfa+l)/2L“(y) 

(23) 

A -tE 

A  P 

(24) 

(25) 

a  =  integer(>l) 

(26) 

y  =  Aexp[-P(r  -  re)J 

(27) 

Nn,aL%  is  the  normalized  associated  Laguerre  polynomial,  r  in  Eq.  27  is  either  R\  or 
rx.  The  parameters  p,  re,  and  De  are  associated  with  the  reduced  mass,  equilibrium 
separation,  fundamental  frequency  and  dissociation  energy  of  the  motion  of  coordinate  r. 
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In  practice  re,  ue,  and  De  are  usually  treated  as  variational  parameters  and  optimized 
accordingly. 

It  is  difficult  to  embed  the  P3  permutation  symmetry  of  ^3-type  molecules  consisting 
of  three  identical  nuclei  in  Tennyson's  method.  However,  the  P2  permutation  symmetry 
of  j4I?2-type  molecules  can  be  easily  built  in45  (also  see  Eq.  18).  Since  the  potential 
energy  function  of  such  a  molecule  is  invariant  under  an  interchange  of  these  two  B  atoms, 
the  Hamiltonian  does  not  couple  the  angular  basis  functions  of  even  j  with  angular  basis 
functions  of  odd  j,  and  we  can  treat  these  two  cases  separately45. 

If  we  treat  an  ^3-type  system  with  only  the  P2  symmetry  embedded  into  the  basis  set 
functions,  the  final  converged  result  should  still  satisfy  the  P3  symmetry.  The  symmetry 
would  manifest  itself  in  the  structure  of  the  eigenenergy  levels  and  in  the  shape  of  the 
eigenfunctions  plotted  in  a  set  of  appropriately  symmetrized  coordinates.  If  an  eigenstate 
obtained  from  the  even  basis  set  is  nondegenerate,  it  must  belong  to  an  Ai  irreducible 
representation  of  P3.  If  an  eigenstate  obtained  from  the  odd  basis  set  is  nondegenerate, 
then  it  generates  an  A2  irreducible  representation  of  P3.  If  one  eigenstate  from  the  even 
basis  set  and  one  eigenstate  from  the  odd  basis  set  are  degenerate  with  each  other,  they 
must  belong  to  an  E  irreducible  representation  of  P3. 

The  Tennyson-Sutcliffe  (TS)  method  has  some  nice  properties.  It  does  not  depend 
on  special  molecular  geometries44,  and  is  applicable  to  motions  with  large  vibrational 
amplitude.  The  choice  of  known  analytic  basis  functions  makes  it  possible  to  obtain  most 
parts  of  the  Hamiltonian  matrix  elements  analytically  and  thereby  save  an  appreciable 
amount  of  computation  time.  The  only  computation-intensive  part  of  the  calculation  is 
the  diagonolization  of  the  Hamiltonian  matrix. 

Numerical  details  and  results 

The  code  we  used  for  the  variational  state  calculation  is  called  TRIATOM,  and  was  ob¬ 
tained  from  the  CPC  Program  Library  of  Queen's  University,  Belfast,  Northern  Ireland45. 
We  initially  made  small  test  runs  on  Sun  workstations  and  a  micro  VAX.  The  major  part 
of  the  calculations  was  done  on  the  SCS-40  mini-supercomputer  of  the  San  Diego  Super¬ 
computer  Center  (SDSC). 

The  ro-vibrational  motion  of  the  ion  has  been  treated  by  Tennyson  et  al.  55-57. 
Since  this  molecule  has  some  resemblance  with  the  Hz  system  in  which  we  are  inter¬ 
ested,  we  repeated  part  of  the  calculation  for  with  total  angular  momentum  J  =  0 
for  gaining  experience  with  this  code.  We  adopted  the  same  values  of  the  parameters 
re,  wr, ,  Dr, ,  /?„,  ujr,  and  Dr,  used  previously55-56,  which  are  listed  in  Table  9.  The 
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Table  9 


Optimized  parameters  of  the  Morse— like  functions  in  R\  and  r*  for  Hf  with  J  =  0. 


De(au) 

ue(au) 

re(au) 

o 

4 

0.230  “ 

0.0085  a 

1.71  ° 

o 

II 

*-> 

0.205  6 

0.0118  b 

2.10  b 

a:  Parameters  for  the  Morse-like  functions  in  R\. 
b:  Parameters  for  the  Morse-like  functions  in  rA. 
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H £  potential  surface  used  in  our  calculation  is  that  included  in  the  TRIATOM  package 
for  code  testing58.  It  is  different  from  the  one  used  in  previous  publications55-57.  As  a 
consequence,  our  results  should  not  be  in  perfect  agreement  with  the  latter.  The  conver¬ 
gence  test  and  the  final  results  are  listed  in  Tables  10  and  11,  which  show  that  the  lowest 
ten  states  are  well  converged,  with  the  size  of  the  largest  basis  set  being  the  same  as  that 
used  previously55-57. 

For  H3,  the  only  currently  available  potential  energy  surface  which  supports  bound 
ro-vibrational  states  of  nuclear  motion  is  the  upper  manifold  of  the  DMBE  surfaces22  for 
the  2 p  2E'  electronic  states,  if  coupling  between  these  two  manifolds  is  neglected.  We 
calculated  the  ro-vibrational  bound  states  with  total  angular  momenta  J  =  0  and  1. 

Initially  we  used  small  sets  of  basis  functions  to  optimize  the  re,  wr>,  Z?r<!,  Re,  ujrc 
and  parameters.  Since  it  had  been  showed  previously  that  the  optimized  parameters 
for  even  j  basis  are  more  or  less  the  same  as  those  for  odd  j  basis45,  we  only  did  the 
optimization  for  the  even  j  basis  calculation.  This  optimization  was  done  with  the  basis 
set  defined  to  be  rnmax  =  8,  nmax  =  8,  jmax  =  16,  and  Ntotai  =  576  for  J  =  0,  and 
mmax  =  6,  nmax  =  6,  jmax  =  15,  and  Ntotai  =  382  for  J  =  1.  For  the  case  of  J  =  0, 
we  put  the  emphasis  on  the  lower  5  states,  while  for  J  —  1  we  considered  the  average 
effect  of  the  parameter  tuning  on  the  lower  12  states.  Because  the  optimization  process 
is  actually  done  manually,  in  a  finite  range  of  the  six  dimensional  parameter  space,  with 
limited  guidance  from  physical  considerations  of  Eqs.  23  to  27,  it  is  possible  that  a  local 
minimum  may  be  accepted  as  the  global  minimum  since  there  is  no  sure  indication  that  the 
global  minimum  has  been  reached.  Fortunately,  the  larger  the  basis  set,  the  less  sensitive 
the  results  are  to  changes  of  those  parameters.  The  results  of  the  optimized  parameters 
are  listed  in  Table  12. 

We  then  increased  the  size  of  the  basis  set  and  tested  the  results  for  convergence. 
We  analyzed  the  importance  of  each  basis  function  for  a  given  basis  size  carefully,  and 
let  the  results  guide  us  to  achieve  an  efficient  way  of  increasing  the  size  of  the  basis  set. 
The  convergence  test  results  for  the  Hz  molecule  are  listed  in  Table  13,  which  shows  that 
the  energy  levels  are  not  well  converged  as  for  .  In  general  the  lower  states  are  better 
converged  than  the  upper  ones.  The  calculation  turned  out  to  be  limited  by  the  amount  of 
computer  memory  we  could  access  at  that  time,  which  was  3  Mwords.  The  final  results  for 
Hz  are  listed  in  Table  14,  with  the  sizes  of  basis  set  in  the  range  of  1100  to  1600  and  SCS-40 
CPU  times  ranging  from  10  minutes  to  30  minutes.  The  results  agree  with  those  from  the 
hyperspherical  coordinate  propagation  method  well,  as  discussed  in  the  next  section. 

One  of  the  most  important  reasons  why  we  need  such  a  large  basis  set  for  the  Hz 
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Table  10 


Convergence  test  of  Hf 


with  J  —  0  and  even  j  basis  functions  for  the  lower  ten  states'*’6. 


#3+. 

with  J  —  0  and  even  j  basis  functions 

m=9,  n=7,  L=  14 

m=10,  n=8,  L=14 

m=ll,  n=8,  L=18 

N  =  340 

N  =  616 

N  =  880 

-7.067955 

-7.067967 

-7.067977 

-6.816079 

-6.816119 

-6.816119 

-6.750378 

-6.750436 

-6.750437 

-6.590602 

-6.590784 

-6.590786 

-6.567921 

-6.568055 

-6.568060 

-6.512097 

-6.512965 

-6.512967 

-6.441278 

-6.441708 

-6.441708 

-6.367264 

-6.367514 

-6.367530 

-6.338421 

-6.338839 

-6.338893 

-6.290215 

-6.290872 

-6.290870 

a  :  The  unit  of  energy  is  104  cm-1. 

b  :  See  reference  45  for  the  definitions  of  (m,  n,  L)  and  N. 
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Table  11:  Hf  J  =  0  bound  state  energies'*. 


Even  basis 

Odd  basis 

Tennyson's6 

1 

Symmetry 

Tennyson's6 

present" 

results 

results 

results 

0.00000 

0.00000 

Ai 

0.24944 

0.25185 

E 

0.24943 

0.25139 

0.31911 

0.31753 

Ai 

0.47250 

0.47718 

Ai 

0.49583 

0.49991 

E 

0.49580 

0.49990 

0.55453 

0.55500 

E 

0.55449 

0.55485 

0.62768 

0.62626 

Ai 

0.69444 

0.70044 

Ai 

0.72350 

0.72907 

E 

0.69433 

0.70030 

Aq 

0.74513 

0.75069 

0.77403 

0.77710 

Ai 

a:  The  energy  Is  in  10*  cm-1  and  its  origin  is  the  ground  state  energy, 
b:  See  references  55-57. 
c:  See  text  and  reference  58. 
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Table  12 


Optimized  parameters  of  the  Morse-like  functions  in  R),  and  rA 
for  with  J  =  0  and  J  —  1. 


De(au) 

w„(au) 

re(au) 

J  =  0 

0.230  a 

0.0130  a 

1.96  “ 

J  =  1 

0.262  a 

0.0100  ° 

2.01  a 

J  =  0 

0.262  6 

0.0122  b 

2.09  6 

J  =  1 

0.232  6 

0.0102  6 

2.32  6 

a:  Parameters  for  the  Morse-like  functions  in  R\. 
b:  Parameters  for  the  Morse-like  functions  in  r*. 
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Table  13 


Convergence  teat  of  H3  with  J  =  0  and  even  j 
basis  functions  for  the  lower  ten  states0,6. 


Hz 

with  J  =  0  and  even  j  basis  functions 

m=15,  n=13,  L=16 

m=16,  n=  13,  L= 18 

10=19,  n=19,  L=26 

N  =  757 

N  =  1067 

N  =  1368 

-0.824614 

-0.826261 

-0.827333 

-0.662336 

-0.664153 

-0.665618 

-0.512776 

-0.514696 

-0.516372 

-0  376072 

-0.377987 

-0.379855 

-0.369457 

-0.369877 

-0.370206 

-0.253486 

-0.256018 

-0.236860 

-0.237305 

-0.134838 

-0.140441 

-0.119736 

-0.120358 

-0.049332 

-0.053411 

a  :  The  unit  of  energy  is  104  cm™1. 

b  :  See  reference  45  for  the  definitions  of  (m,  n,  L)  and  N. 
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Table  14 


Bound  state  energies  without  consideration  of  the  geometric  phase". 


J  =  06 

J  = 

L  even  parity*" 

J  = 

1  odd  parityc 

3.7210 

Ai 

3.7218 

3.7282 

A2 

3.7294 

3.7264 

E 

3.7276 

3.9216 

Ax 

3.9223 

3.9284 

A2 

3.9297 

3.9266 

E 

3.9281 

4.1067 

Ax 

4.1073 

4.1130 

A2 

4.1145 

4.1114 

E 

4.1131 

4.2759 

Ai 

4.2766 

4.2817 

a2 

4.2839 

4.2802 

E 

4.2831 

4.4282 

Ax 

4.4301 

4.4336 

A2 

4.4386 

4.4322 

E 

4.4398 

4.5621 

Ax 

4.5734 

4.5665 

A2 

4.5803 

4.5656 

E 

4.5894 

4.2886 

E 

4.2886 

4.2955 

E 

4.2956 

4.2971 

Ax 

4.2975 

4.2969 

A2 

4.2972 

4.2904 

E 

4.2908 

4.4533 

E 

4.4533 

4.4596 

E 

4.4598 

4.4610 

Ax 

4.4618 

4.4608 

A2 

4.4615 

4.4550 

E 

4.4557 

4.5980 

E 

4.5983 

4.6036 

E 

4.6048 

4.6049 

Ax 

4.6083 

4.6047 

a2 

4.6093 

4.5996 

E 

4.6028 

4.7212 

E 

4.7261 

E 

4.7349 

4.7272 

Ax 

4.7370 

4.7270 

a2 

4.7355 

4.7225 

E 

4.6806 

Ax 

4.6813 

4.6871 

A2 

4.6893 

4.6842 

E 

4.6878 

a:  The  energy  is  in  eV  and  its  origin  corresponds  to  the  bottom  of  the  ground  electronic  state  of  the 
isolated  H2  molecule. 

b:  The  left  column  gives  the  results  of  the  hyperspherical  coordinate  propagation  method  and  the  right  col¬ 
umn  the  TS  method  results.  The  central  column  gives  the  irreducible  representation  of  the  permutation 
group  of  the  nuclei  to  which  the  spacial  part  of  the  nuclear  wavefunction  belongs. 


system  is  that  the  upper  manifold  of  the  Hz  2 p  2E'  surface  has  a  cone-shape  structure  with 
a  sharp  tip  at  the  equilateral  triangle  nuclear  geometry,  while  the  surface  has  a  nice 
smooth  round  shape  instead.  The  Morse-type  basis  set  functions  are  good  for  representing 
the  harmonic-like  oscillation  of  the  ion.  However  they  are  not  good  for  representing 
the  motion  on  such  a  conicaily  shaped  surface.  As  a  result,  more  basis  functions  are  needed 
to  obtain  converged  results  i/3  than  for  ,  if  the  same  degree  of  convergence  is  required. 
For  the  2s  2A[  and  2 px  2A'{  electronic  states,  the  /J3  is  really  like  a  core  plus  a  Rydberg 
electron  in  2s  and  2 p2  states  respectively,  which  interact  weakly  with  the  Hf  core.  For 
this  reason  the  snapes  of  the  potential  energy  surfaces  of  those  two  electronic  states  should 
be  simiKr  to  that  of  the  Hf  in  its  ground  electronic  state.  This  suggests  that  a  small  basis 
sets  (01  size  about  800)  should  give  converged  results  for  the  2s  2A\  and  2 pz  2A'^  potential 
energy  surfaces. 

Finally,  let  us  consider  the  shapes  of  the  ro-vibrational  wavefunctions  to  see  if  the 
final  converged  calculations  yield  wavefunctions  with  the  right  P3  symmetry.  We  ploted 
the  wavefunctions  in  a  system  of  symmetrized  coordinates59.  Figure  7  contains  contour 
lines  of  the  wavefunctions  for  Hf  with  total  angular  momentum  J  —  0  and  basis  set  size 
N=880,  and  Figure  8  for  Hz  with  J  =  0  and  basis  set  size  N=1363.  The  plots  show  that  the 
wavefunctions  do  not  display  P3  (».«.,  Czv )  symmetry,  even  for  the  highly  converged  state 
of  .  The  reason  for  this  behavor  is  that  in  general,  the  convergence  of  the  eigenvector 
is  poorer  than  that  of  the  eigenvalue  in  a  numerical  eigenvalue-eigenvector  problem.  In 
order  to  get  the  right  symmetry  with  a  reasonable  basis  set  size,  the  symmetry  has  to 
be  embeded  into  this  basis  set  before  the  variational  calculation  is  performed.  This  is 
difficult  to  do  with  Tennyson’s  code,  so  we  developed  a  new  method  to  achive  this,  which 
is  described  in  section  3.2. 
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Figure  7:  Contour  plot  of  the  wavefunction  for  the  lowest  J  =  0,  A^-type 
ro-vibrational  state,  in  symmetrized  hyperspherical  coordinates50.  Depicted  is  a  cut  at 
constant  Yx  chosen  to  be  2.1  bohr,  for  which  the  potential  energy  function  of  Hf  has 
a  minimum  (at  Xx  =  Zx  =0).  The  maximum  (near  the  center  of  the  plot)  of  the 
wavefunction  was  set  equal  to  1.0,  and  contours  for  #  =  0.9  to  0.1  in  steps  of  0.1  are 
displayed.  The  Xx  and  Zx  coordinates  are  in  bohr. 


Figure  8:  Contour  plot  of  the  wavefunction  'P  for  the  lowest  J  =  0,  Ai-type  H3 
ro-vibrational  state,  in  symmetrized  hyperspherical  coordinates59.  Depicted  is  a  cut  at 
constant  Y\  chosen  to  be  1.967  bohr,  for  which  the  potential  energy  function  of  H3  has 
a  minimum  (at  X\  =  Z\  =0).  The  maximum  (near  the  center  of  the  plot)  of  the 
wavefunction  was  set  equal  to  1.0,  and  contours  for  =  0.9  to  0.1  in  steps  of  0.1  are 
displayed.  The  X\  and  Z\  coordinates  are  in  bohr. 
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3.2  Hyperspherical  Coordinate  Propagation  Method 


The  hypersphericai  coordinate  approach  has  been  successfully  used  in  recent  years 
for  the  calculation  of  bound  ro-vibrational60-61  and  scattering62-07  states  of  triatomic 
systems.  Our  motivations  of  using  it  to  study  the  bound  ro-vibrational  motion  of  H3  in 
its  upper  manifold  of  the  electronic  2p  2E'  states  are  three-fold: 

1.  The  full  P3  nuclear  permutation  symmetry  of  the  identical  triatomic  system  can  be 
easily  embedded  into  the  basis  set  functions  so  that  the  final  wavefunctions  of  the 
ro-vibrational  states  will  have  the  correct  symmetry.  This,  as  seen  in  the  previous 
section,  is  difficult  to  achieve  with  the  Tennyson-Sutcliffe  variational  method. 

2.  The  ground  electronic  state  2 p  2£'(l)  of  H3  has  a  conical  intersection  with  the  first 
excited  electronic  state  2p  2E'{ 2)  at  equilateral  triangle  nuclear  configurations19,22. 
The  Born-Oppenheimer  real  electronic  wavefunctions  undergo  a  sign  change  when  one 
follows  a  close  path  in  nuclear  configuration  space  around  the  line  along  which  the 
two  states  conically  intersect 19,22,49-52 .  Since  the  total  electro-nuclear  wavefunction 
is  continuous  and  single  valued,  there  has  to  be  a  compensating  sign  change  on  the 
nuclear  part  of  the  wavefunction49-52.  This  is  referred  to  in  the  literatures  as  the 
Molecular  Aharonov-Bohm  (MAB)  effect51-52.  This  effect  will  modify  the  energy 
levels  of  the  bound  ro-vibrational  states  of  the  upper  electronic  state  2 p  2E[ 2)  which 
would  exist  in  the  absence  of  interaction  with  the  ground  state.  It  is  easy  to  take  this 
MAB  effect  into  account  with  the  symmetrized  hyperspherical  coordinate  propagation 
method  described  below,  while  this  is  not  possible  with  the  TS  variational  method 
in  its  present  format.  Our  results  should  help  to  solve  the  controversy  as  to  whether 
these  ro-vibrational  bound  states  have  been  experimentally  detected  or  not2-4. 

3.  The  corresponding  wavefunctions  will  be  needed  to  calculate  the  overlap  integrals 
with  the  ground  state  scattering  wavefunctions  (see  Eqs.  1  and  2).  Since  we  calculate 
the  latter  using  hyperspherical  coordinates,  it  makes  sense  to  obtain  the  bound  state 
wavefunctions  with  the  same  coordinates. 

In  the  absence  of  coupling  between  these  two  electronic  states,  the  excited  2 p  2E'( 2) 
surface  should  support  bound  ro-vibrational  states.  Initially,  without  including  the  MBA 
effect,  we  computed  those  states  using  a  symmetrized  hyperspherical  coordinate  propa¬ 
gation  method,  and  compared  the  results  with  those  of  the  variational  method  (TS),  for 
total  angular  momentums  7  =  0  and  7  =  1.  Then  we  did  the  calculation  including 
the  MAB  effect,  and  obtained  completely  different  energy  levels68.  The  MAB  effect  also 
modifies  the  scattering  calculation  in  H  +  H2  system69. 
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The  Method 


Let  Aa,  Ap ,  A-,  be  the  atoms  of  the  system,  and  (A,i/, /c)  be  any  cyclic  permutation 
of  (a,/?, 7).  T\  is  the  mass-scaled70  internuclear  vector  for  the  diatom  AUAK  and  R>  the 
mass-scaled  vector  of  A\  with  respect  to  the  center  of  mass  of  A^A*.  The  hyperspherical 
method  uses  the  hyper-radius  p  =  (Rf  +  r2)  s  to  describe  the  global  size  of  the  triatomic 
system  and  a  set  of  five  angles  f  to  describe  its  shape  and  orientation  in  space62-67,70.  In 
the  Born-Oppenheimer  approximation,  the  electro-nuclear  wavefunction  can  be  written  as 
a  product  of  the  electronic  part  ipe,  which  we  choose  to  be  real,  and  the  nuclear  part.  The 
latter  can  be  factored  into  a  nuclear  spin  part  and  a  spatial  part  r/;JMnr .  J  is  the  total 
nuclear  angular  momentum  quantum  number,  M  its  projection  onto  a  laboratory-fixed  z 
axis,  Ft  the  parity  with  respect  to  inversion  of  the  nuclear  coordinates  through  their  center 
of  mass  and  F  the  irreducible  representation  of  the  nuclear  permutation  group  (P3)  to 
which  q/JMnr}  the  electro-nuclear  wavefunction  excluding  the  nuclear  spin  part,  belongs: 

qJMUT  _  0JA/nr(p,<r)0e(ge;p,f)  (28) 

qe  refers  to  the  set  of  all  electronic  coordinates  (spacial  and  spin).  xl>JMXir  is  an  eigenfunc¬ 
tion  of  the  nuclear  motion  Hamiltonian 


k*  sa  .a  a2 

2ft  dp  dp  2  up2 


+  V(p,f) 


(29) 


where  p  is  the  three-body  reduced  mass,  A  the  grand  canonical  angular  momentum  and  V 
the  Bom-Oppenheimer  electronic  potential  energy  function.  The  nuclear  function 
is  expanded  in  a  basis  of  local  hyperspherical  surface  functions  (LHSF)  $^wnr  ; 

v>jMnr(p,o  =  4  (30) 

O  3  - 

r  n 

The  LHSF  are  defined  as  the  eigenfunctions  of  the  nuclear  Hamiltonian  at  fixed  hyperradius 
P ■ 

15^5  +  (31) 

The  coefficients  F/nr  in  Eq.  30  are  solutions  of  a  set  of  coupled  differential  equations- 
in  p,  whicl  ve  solve  using  piece-wise  diabatic  bases62-67.  For  assumed  values  of  the 
rovibrational  energies,  the  solutions  are  propagated  forward  and  backward  from  small  and 
large  p  values  where  they  have  negligible  amplitudes.  The  energy  is  scanned  iteratively 
until  the  forward  and  backward  solutions  match  smoothly  at  an  intermediate  value  of  p. 
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In  the  present  study,  we  use  the  Whitten-Smith71-72  definition  of  the  five  angular 
coordinates  f  as  modified  by  Johnson73.  Three  Euler  angles  ( a(3^ )  specify  the  orientation 
of  the  body  frame  in  space.  The  axes  of  this  frame  lie  along  the  principle  axes  of  inertia. 
The  Z  axis  is  parallel  to  T\  x  R>  and  the  X  axis  is  associated  to  the  smallest  moment  of 
inertia  and  is  oriented  such  that  rXx  >  0.  Two  angles  (0,<px)  describe  the  shape  of  the 
molecular  triangle  and  are  defined  by: 


7T  0 ,  .  ,<P\, 
rxx  =  pcos(-  -  -)sin(— ) 

(32) 

r\y  =  — psin(—  -  -)cos(— ) 

(33) 

Rxx  =  pcos(^  -  ^)cos(-y ) 

(34) 

Rxy  =  Psin(~  -  ^)sin(-) 

(35) 

The  ranges  for  these  angles  are  0  <  0  <  \  and  0  <  <px  <  2n r.  0  =  0  corresponds  to  the 
symmetric  top  configuration  (an  equilateral  triangle  for  three  identical  particles)  in  which 
the  principal  axes  of  inertia  X  and  Y  are  undefined. 

The  grand  canonical  angular  momentum  is  given  explicitly  by67-69: 


P  =  -4A2{- 


9  .  a 

.  sm20—  + 

sm20  d0  dO 


+  fcSl  + 


n 


cos2  0  sin20 


sin20  d<pl 

+  +  p] 

+  cos*»l  ++  -1 


d 2  ,  4ihcos0  ?  d 

}  +  —rrrJ*T— 

sin20  oipx 


(36) 


A  A  A 

where  Jz  is  the  body-fixed  Z  component  of  the  total  angular  momentum  J,  and  J±  = 
Jx  ±  iJy  • 

Eq.  31  is  solved  variationally  by  expansion  in  a  body-fixed  basis  Xn^n^  built  with 
products  of  simple  analytical  functions60  : 


D3mk  is  a  Wigner  rotation  matrix53  and  nv  is  integer  or  half  of  an  odd  integer.  fne  (0) 
are  simple  trigonometric  functions,  such  that  the  LHSF  have  correct  behaviors  near  the 
singularities  of  the  kinetic  energy  operator  0  =  0  and  In  practice,  the  /n#  can  be  chosen 
as  the  functions  cos (n#0)  or  sin(ri00),  with  no  integer  or  half  odd  integer,  in  terms  of  which 
the  hyperspherical  harmonics  (whose  0  dependence  is  usually  written  as  a  polynomial  in 
cos(0))  can  be  written. 
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We  now  focus  attention  on  the  special  case  of  three  identical  nuclei  and  describe 
how  to  build  electro- nuclear  wavefunctions  \^J’wnr  which  are  bases  for  the  irreducible 
representations  of  the  permutation  group  of  the  nuclei  (P3).  The  operations  of  this  group 
correspond  to  simple  changes  in  <p\  (which  are  related  to  the  isomorphism  between  P3  and 
C3„)  as  indicated  in  Table  15.  If  ±1)  is  the  symmetry  of  the  electronic  wavefunction 

with  respect  to  the  1/  <-►  x  permutation,  then  the  linear  combinations  defined  by 

JMKtl*  JMK  ,  en  e  (_-,\J+K+2n„  JM,-K 

give  electro-nuclear  wavefunctions  ijrjwnr  (Eq.  28)  having  e*”(=  ±1)  symmetry  with 
respect  to  the  u  k  permutation. 

If  there  is  no  conical  intersection  between  electronic  states,  the  electronic  wavefunc¬ 
tion  Tl)e{qe\ P,$)  belongs  to  a  one  dimensional  representation  of  the  nuclear  permutation 
group  (Ai  for  =  +1,  or  A 2  for  =  —1).  Table  16  indicates  how  the  total  angular 
momentum,  the  parity  and  the  irreducible  representation  T  of  P3  to  which  \J/JMnr  belongs 
determines  the  set  of  quantum  numbers  nv. 

The  hyperspherical  method  uses  20  n0  values,  between  4  (  Ai  or  A2  symmetry  )  and 
8  (  E  symmetry  )  \nv\  values  (Eqs.  37  and  38),  and  between  6  (  At  or  A2  symmetry 
)  and  12  (  E  symmetry  )  LHSF  (Eqs.  30  and  31).  The  LHSF  have  been  computed  at 
typically  50  p  values  between  1.5  bohr  and  6.5  bohr.  The  convergence  of  the  LHSF  and 
rovibrational  energies  is  of  the  order  of  10-4  eV.  The  compactness  of  the  hyperspherical 
expansion  comes  from  the  fact  that  the  potential  energy  around  the  Y\  axis  ( 0  =  0)  is 
nearly  cylindrical  (small  number  of  n\v\  values)  and  from  the  steep  increase  of  the  the 
potential  as  a  function  of  6  (small  number  of  LHSF). 

Results  and  Conclusions 

Figure  9  illustrates  the  main  features  of  the  electronic  potential  energy  surface  in 
an  internal  configuration  space50.  The  energy  levels  obtained  from  the  calculation  in  the 
hyperspherical  coordinates  are  converged  to  within  1  cm-1,  and  are  very  close  to  those 
gotten  by  using  the  TS  variational  method,  with  the  former  one  more  converged  and  more 
accurate.  The  hyperspherical  coordinates  propagation  method  makes  the  symmetry  as¬ 
signment  of  these  levels  much  easier.  Table  14  compares  the  ro-vibrational  energy  levels 
obtained  by  these  two  methods  excluding  the  MAB  effect.  Table  17  gives  the  energy  levels 
when  the  MAB  effect  is  included,  obtained  from  the  hyperspherical  coordinate  propaga¬ 
tion  method  only.  Comparison  between  Tables  14  and  17  shows  that  this  effect  is  very 
important.  Figure  10  is  the  graphic  representation  of  Tables  14  and  17.  (/xt,  p2,  0  are 
the  quantum  numbers  assigned  to  those  states68. 
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Table  15 


Effect  of  permutations  of  the  nuclei  on  the  angle  <px. 


Permutation 

Pa„/ 

P.«xb  ■ 

PnXuC 

P  d 

*  UK 

Pxud 

Pa/ 

Value  of  <pxe 

<fix 

'PA  +  ^f 

2t  -  <f>x 

'f-Vx 

a  :  Px uk  i3  the  identity  permutation. 

b  :  PUk\  refers  to  the  cyclic  permutation  Xi/k  — ►  ukX. 
c  :  Pk\u  refers  to  the  cyclic  permutation  Xi/k,  — *  kXv. 
d  :  Pt]  refers  to  the  pairwise  permutation  of  nuclei  t  and  j. 

e  :  The  changes  in  <px  are  true  modulo  2ir,  since  <px  must  remain  in  the  range  [0,27rJ. 
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Figure  9:  The  two  dimensional  plot  of  the  upper  sheet  of  the  DMBE  potential 
surfaces22.  The  hyperspherical  coordinates  are  used59.  The  Zx  value  is  zero.  The  conical 
intersection  of  those  two  electronic  states  happens  along  the  Y\  axis.  The  equipotentials 
are  equally  spaced  by  0.25  eV  from  3.0  eV  to  5.0  eV.  Xx  and  Y\  are  given  in  bohr. 


Table  16 


Choice  of  nv  for  each  parity  n  and  irreducible  representation  T 
of  the  nuclear  permutation  group  P3. 


n 

rc 

Even, without  phase3 

Ax/A2 

Odd, with  phase*1 

E 

Even, with  phaseb 

A1/A2 

Odd, without  phase3 

E 

3  md 
3m±ld 


3m  + 

3  m±±a 


a  :  without  consideration  of  the  geometric  phase  due  to  the  conical  intersection, 
b  :  with  consideration  of  the  geometric  phase  due  to  the  conical  intersection, 
c  :  T  is  the  irreducible  representation  of  P3  to  which  yJMnr  belongs, 
d  :  m  is  a  non-negative  integer. 
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Figure  10:  Ro-vibronic  energy  levels  associated  to  the  upper  sheet  of  the  DMBE 
potential  surfaces  of  H322.  The  full  lines  are  the  levels  including  the  effect  of  the  geometric 
phase  while  the  dashed  ones  exclude  that  effect.  The  quantum  numbers  vlt  v2  and  l  are 
defined  in  reference  68.  The  origin  for  the  energy  scale  is  the  bottom  of  the  isolated  ground 
electronic  H 2  potential  energy  curve.  These  levels  are  for  the  J  =  0  states,  but  the  J  =  1 
levels  are  nearly  degenerate  with  them,  the  splitting  being  of  the  order  of  10-2  eV.  Their 
nuclear  permutation  symmetries  depend  on  J  and  on  the  parity  IT,  as  well  as  whether  the 
geometric  phase  is  or  is  not  included  (see  Table  15  and  Table  16).  There  are  two  levels  for 
each  of  the  sets  of  quantum  numbers  (t>i  =  0,v2  =  /  =  |)  and  (t»i  =  l,v2  =  /  =  |),  which 
would  be  degenerate  if  the  potential  were  exactly  cylindrically  symmetric  around  the  Y\ 
axis  (see  text  and  Figure  9). 
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We  have  described  in  this  section  a  new  hyperspherical  coordinate  propagation  method 
for  the  calculation  of  bound  ro-vibrational  states  of  triatomic  system.  This  method  is  well 
adapted  to  systems  of  three  identical  particles,  because  it  allows  easy  inclusion  of  the  full 
permutation  symmetry  of  the  system  and  of  the  effect  of  conical  intersections  on  the  phase 
of  the  nuclear  wavefunction  (MAB  effect).  Since  there  is  not  MAB  effect  for  the  2s  2A\ 
and  2p3  electronic  potential  energy  surfaces,  we  can  use  both  methods  in  the  study 
of  their  bound  ro-vibrational  states.  The  TS  method  can  easily  furnish  the  approximate 
locations  of  the  ro-vibrational  eigen-energies,  while  the  hypersperical  one  permits  us  to 
scan  energies  near  those  approximate  ones  and  obtain  quickly  more  accurate  levels  and 
their  corresponding  wavefunctions,  having  the  exact  P3  symmetry  built  in.  This  work  has 
beer,  published  recently68. 
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4.  Scattering  Wavefunction 


4.1  General  Considerations 

In  the  theoretical  treatment  of  predissociation  processes11,12  and  of  transition  state 
spectroscopy  (TTS)39,  the  unbound  scattering  wavefunction  with  a  proper  physical  bound¬ 
ary  condition  and  normalization  is  required. 

In  a  typical  quantum  scattering  calculation,  the  scattering  matrix  S  is  obtained  instead 
of  the  physical  wavefunction  with  proper  boundary  conditions  and  normalization.  The 
reason  is  that  the  scattering  matrix  S  contains  all  the  information  about  the  state-to- 
state  scattering  cross  sections  (which  can  be  compared  with  experimental  results),  and 
therefore  makes  the  construction  of  the  physical  wavefunction  unnecessary.  Many  methods 
of  scattering  calculation  take  advantage  of  this  fact  and  obtain  the  S  matrix  without 
constructing  the  physical  wavefunctions62-67,69.  Only  a  few  calculations  involved  with 
scattering  wavefunctions  have  been  published  for  triatomic  systems,  both  for  collinear74 
and  3D75  cases. 

Because  the  S  matrix  involves  the  ratio  between  the  the  in-coming  parts  and  the  out¬ 
going  ones  of  the  scattering  wavefunctions  in  the  asymptotic  regions  (where  the  triatomic 
system  becomes  a  free  atom  plus  a  diatomic  molecule,  and  the  wavefunctions  can  be  ex¬ 
pressed  as  a  sum  of  products  of  the  ro-vibration  wavefunctions  of  the  diatomic  molecule 
and  the  plan-wave  functions  of  the  free  atom  with  respect  to  the  center  of  mass  of  the 
diatomic  molecule) ,  it  allows  the  determination  of  the  physical  wavefunction  in  the  asymp¬ 
totic  regions  for  a  given  choice  of  the  in-coming  parts.  This  furnishes  the  explicit  boundary 
problem  of  the  Schrodinger  equation.  The  solution  will  correspond  to  the  real  physical  pro¬ 
cess  of  interest 

Using  the  coupled-channel  method  in  hyperspherical  coordinates,  the  wavefunction 
can  be  expressed  as  am  expansion  in  basis  set  functions  of  the  hyperangles  with  coefficients 
as  functions  of  the  hyperradius,  the  Schrodinger  (partial  differential)  equation  can  be 
transformed  into  a  set  of  second  order  ordinary  coupled  equations  as 

(1^  +  Q(x)]¥(x)  =  0  (39) 

with 

Q(x)  =  (2M/ft2)[£I-V(x)]  (40) 

where  I  is  the  identity  matrix,  and  V(x)  the  interaction  matrix.  E  is  the  known  total  energy 
of  the  system,  't'(x)  the  expansion  coefficient  matrix62-67,76.  Here  x  is  the  hyperradius. 
Detail  of  the  definitions  of  those  matrices  cam  be  found  elsewhere62-67,76 
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Once  scattering  matrix  S  has  been  determined  by  a  forward  integration  method,  the 
physical  wavefunction  'l'(x)  can  be  calculated  from  a  backward  integration  of  Eq.  39.  The 
initial  conditions  are  obtained  by  projecting  the  wavefunction  in  the  a®-  ptotic  region, 
which  is  obtained  from  S,  into  the  hyperspherical  surface  function  at  large  hyperradius. 
The  backward  integration  is  then  conducted  to  a  small  hyperradius  where  the  interaction 
is  strongly  repulsive  and  the  wavefunction  vanishes. 

Previous  calculations62-64  are  all  based  on  the  logarithmic  derivative  method77  which 
is  not  suitable  for  the  construction  of  the  physical  wavefunction  in  the  backward  propaga¬ 
tion  because  it  furnishes  the  logarithmic  derivative  matrix  rather  than  the  wavefunction 
matrix  itself.  We  have  chosen  the  renormalized  Numerov  method78  for  that  backward 
integration. 

4.2  Renormalized  Numerov  Method 

In  the  renormalized  Numerov  method,  the  Eq.  39  is  transformed  into 

[I  -  Tn+1]$n+1  -  [21  -  Tn]*n  +  [I  -  =  0  (41) 

by  using  finite  difference  scheme,  where 

*n  =  ¥(*»)  (42) 

and 

T,  =  -(/i2/l2)Q(x„)  (43) 

Here  h  is  the  spacing  between  the  N  +  1  equally  spaced  grid  point  (xo,  x*,  . .  ,  xjv)  and 
Q(x)  is  defined  by  Eq.  40.  The  Fn  matrix  is  defined  as 

F»  =  [I  —  T„]^„  (44) 

If  we  substitute  Eq.  44  back  into  Eq.  41,  we  have 

Fn+i  -  U„Fn  -I-  F„_i  =  0  (45) 

with 

U„  =  (I-T*)-1  (2I  +  10T„)  (46) 

Next,  a  ratio  matrix  R  is  defined  as 

Rn  =  Fn+iF-1  (47) 
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I 

I 

This  transforms  Eq.  45  into  the  two  term  recurrence  relation 

j  Rn  =  U„  -  R"^  (48) 

i 

This  is  the  basic  propagation  relation  for  this  method. 

At  large  hyperradii  (say,  xo  and  Xi),  the  wavefunction  is  known,  so  are  the  expansion 
coefficient  matrices  'I'o  and  which  are  obtained  by  the  wavefunction  projection  in  the 
asymptotic  regions.  The  matrices  Fo  and  Fi  are  obtained  by  using  Eq.  44,  and  then  Ro 
by  using  Eq.  47.  Then  the  backward  propagation  can  be  started  following  Eq.  48,  and  all 
be  generated  accordingly. 

We  have  implemented  the  renormolized  Numerov  propagation  scheme  in  the  collinear 
if3  scattering  calculation  on  its  ground  potential  energy  surface  in  order  to  test  our  version 
of  this  propagator.  The  cc  le  is  currently  being  debugged  and  compared  with  previous 
results  for  forward  integration'9.  Once  this  Ls  completed,  the  backward  integration  method 
for  d  termined  the  physical  wave  function  will  be  extended  to  three  dimensions. 
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5.  Decay  Mechanisms  for  the  //3  System 

The  metastable  //3  molecule  in  the  2 pz  2A2  electronic  state  can  decay  via  radiative 
processes  to  the  lower  2s  2A\,  2 p  2E'{ 2)  and  2p  2E'{\)  electronic  states  and  by  the  predis¬ 
sociation  to  the  unbound  2 p  2E'{  1)  ground  electronic  state.  In  this  section,  we  discuss  the 
calculation  of  the  lifetime  via  each  of  these  decay  channels. 

5.1  Radiative  Decay  and  Selection  Rules 

By  coupling  to  the  radiation  field,  the  H$  molecule  can  emit  radiation  as  it  undergoes 
a  transition  from  an  upper  state  to  a  lower  one.  The  general  theory  of  the  interaction 
between  a  molecular  system  and  a  radiation  field  is  beyond  the  scope  of  this  report  and  can 
be  found  elsewhere17.  We  will  consider  here  three  important  kinds  of  transition:  electric 
dipole  moment  (El),  magnetic  dipole  moment  (Ml)  and  electric  quadrupole  moment 
(E2)  The  lifetime  r  associated  to  each  is  inversely  proportional  to  the  square  of  the 
corresponding  transition  moments: 


Tx  oc  |  |  X  |  *')  |  2  (49) 

where  'F‘  and  ty*  are  the  initial  and  final  wavefunctions  for  the  molecule  and  X  the 
transition  moment  operator  under  consideration  (which  can  be  El,  Ml  or  E2). 

The  equilibrium  geometries  of  the  bound  ro-vibrational  states  of  the  2s  2A\,  2pz  2  A" 
and  the  upper  manifold  of  the  2p  2E'  states  (in  the  absence  of  John- Teller  distortions19,27) 
are  of  that  of  an  equilateral  triangle.  For  this  reason,  the  wavefunctions  of  those  electronic 
states  display  (to  zeroth  order)  symmetry.  Using  this  symmetry,  and  with  the  assump¬ 
tion  that  the  transition  dipole  moments  do  not  change  with  nuclear  geometry,  Herzberg 
d.al.  have  discussed  the  selection  rules  for  the  electric  dipole  radiative  transitions  among 
those  four  electronic  states5,7.  We  have  extended  their  considerations  of  the  electric  dipole 
transition  (El),  to  that  of  the  magnetic  dipole  (Ml),  and  electric  quadrupole  (E2). 

We  use  the  same  assumption  that  the  multipole  moments  do  not  change  with  nuclear 
geometry.  Under  this  condition,  the  integral  in  Eq.  49  can  be  separated  into  a  electronic 
part  and  a  nuclear  part,  which  permits  us  to  rewrite  Eq.  49  as: 

rx  a  |  {%  |  X  |  *£)  f2|  (K  |X|  #i>  f  *  (50) 

where  'F(.  and  stand  for  the  electronic  and  nuclear  wavefunctions  respectively.  The 
electric  dipole  moment  El  (being  a  vector)  forms  an  A 2  +  E'  representation  of  the  D^h 
group,  the  magnetic  dipole  moment  Ml  (being  a  skew  vector)  forms  an  A'2  +  E"  one, 
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and  the  electric  quadrupole  moment  E2  (as  a  second  rank  tensor)  an  A\  +  E1  +  E"  one. 
The  symmetry  requires  that  for  the  integral  over  electronic  coordinates  in  the  right  hand 
side  of  Eq.  50  not  to  vanish,  it  is  necessary  that  the  integrand  (the  product  of  the  initial 
electronic  state,  the  moment  operator  and  the  final  electronic  state)  in  that  integral  must 
have  an  A\  component.  Using  the  properties  of  the  D^h  group,  this  furnishes  the  selection 
rules  for  the  transitions  between  those  four  electronic  states. 

For  the  electric  transition  dipole  moment  (El),  the  transition  can  have  non-zero  intensity 

•  between  the  2 pz  2A%  and  2s  2A\  states, 

•  and  between  the  2s  2A\  and  2 pIiS/  2E'[  1,2)  states. 

For  the  magnetic  transition  dipole  moment  (Ml),  it  can  have  non-zero  intensity 

•  between  the  2 p3  2A 4'  and  2 px%y  2J5#(1,2)  states. 

For  the  electric  transition  quadrupole  moment  (E2),  it  can  have  non-zero  intensity 

•  between  the  2s  2A\  and  2 px,v  2E'[  1,2)  states, 

•  and  between  the  2p2  2A'{  and  2pz,„  2E'{  1,2)  states. 

Since  we  have  calculated  the  potential  energies  of  the  2 pz  2Ar2  and  2s  2A\  states  and 
the  electric  transition  dipole  moment  between  them,  we  can  estimate  the  lifetime  of  this 
decay  process.  Our  estimation  was  done  at  R  =  1.64  bohr  (iZ  being  the  length  of  the  side 
of  the  equilateral  H3)  where  the  nuclear  configuration  is  very  close  to  the  minima  of  both 
potential  energy  surfaces,  using 

=  2.02586  x  10"6(A£)3T?/  (51) 

where  A E  is  the  transition  energy  (A E  =  Ei  —  Ej,  in  cm-1),  and  T*,/  is  the  transi¬ 
tion  dipole  moment  (in  atomic  unit)  between  the  initial  state  (2s  2A\)  to  the  final  state 
(2 pz  2A‘2),  and  the  dipole  transition  lifetime  r*,/  is  in  seconds.  Our  calculated  transition 
energy  is  1299  cm-1  and  the  corresponding  electric  transition  dipole  moment  is  2.68  a.u. 
which  yields  lifetime  of  about  31  ps.  This  compares  with  the  results  of  1422  cm-1,  2.69 
a.u.  and  23  ps  from  Petsalakis  et  al.32  ,  and  1988  cm-1,  2.85  a.u.  and  7.7  ps  from  King 
and  Morokuma27.  The  lifetime  is  very  sensitive  to  the  transition  energy.  If  the  experimen¬ 
tal  transition  energy  (of  993  cm-1)  is  used  (which  includes  the  difference  between  the  zero 
point  energies  of  the  nuclei  motion  on  the  2s  2A[  and  2 pz  2A'2  potential  energy  surfaces), 
the  resulting  lifetime  estimations  are  70  ps  (this  calculation),  62  ps  (  King  and  Morokuma) 
and  70  ps  (Petsalakis  et  al.)  and  all  agree  with  the  experimental  estimated  value  reported 
as  being  in  excess  of  40  ps14. 
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When  the  ro-vibrational  motion  of  the  nuclei  distorts  the  shape  of  the  molecule  away 
from  equilateral  triangle,  the  D^u  is  no  longer  the  symmetry  group  of  the  electronic  wave- 
function,  so  that  those  selection  rules  should  not  be  observed  exactly.  For  a  molecular 
system,  the  strength  of  the  coupling  with  the  radiation  field  generally  decreases  in  the 
order  l)  electric  transition  dipole,  2)  magnetic  transition  dipole,  3)  electric  transition 
quadrupole,  and  so  on.  Only  if  some  selection  rule  prevents  the  stronger  coupling  to  be 
non-zero,  can  the  weaker  coupling  have  a  chance  to  show  its  contribution.  When  Hz 
does  not  have  equilateral  triangular  geometry,  the  electric  transition  dipole  coupling  can 
contribute  to  each  of  the  possible  transition  between  those  lowest  four  electronic  states. 
Therefore,  it  is  appropriate  for  us  to  study  the  electric  dipole  transition  first.  With  the 
potential  energy  surfaces  and  electric  transition  dipole  moments  obtained  by  the  electronic 
calculation  described  in  section  2,  the  nuclear  ro-vibrational  bound  states  can  be  obtained 
by  the  methods  described  in  section  3,  and  then  the  radiative  lifetime  of  2 pz  2 A — +  2s  2A\ 
can  be  calculated. 


5.2  Predissociation 

Predissociation  can  occur  as  a  result  of  coupling  between  two  Bom-Oppenheimer 
electronically  adiabatic  potential  energy  surfaces  associated  with  two  electronic  states  of  a 
molecular  system.  The  Bom-Oppenheimer  wavefunctions  of  the  system  are  not  in  this  case 
true  eigenstates  of  the  total  molecular  Hamiltonian.  Interaction  between  the  quasi-bound 
states  of  the  nuclear  motion  of  an  upper  surface  and  the  unbound(».e.,  scattering)  ones  of 
a  lower  surface  can  cause  a  quantum  leakage  from  the  former  to  the  latter15-16 

5.2.1  Fano's  Theory16 

In  the  simplest  case,  there  are  two  quantum  states  associated  with  the  Hamiltonian 
H.  One  is  bound  and  denoted  as  J  <f>n)  and  the  other  one  unbound  and  denoted  as  |  E). 
They  are  not  exact  eigenstates  of  H  and  satisfy  the  following  conditions: 


(4>r, 


(<t>. 


(<£n  |  <t>n)  =  1 

(52) 

(4>n  |  E)  =  0 

(53) 

(E  |  El)  =  6{E  -  E') 

(54) 

1  H  |  <t>n)  =  En 

(55) 

\H\E')  =  ES{E  -  E') 

(56) 

,\H\E)=  Vn(E) 

(57) 

where  n  designates  the  set  of  quantum  number  which  label  the  bound  state,  and  E  is  the 
energy  for  the  unbound  state.  Vn(E)  is  the  coupling  between  the  bound  state  and  the 
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unbound  state  and  is  usually  very  small.  We  expand  the  true  eigenstate  of  H  as 


|  *„(£)>  =  An(E)  \<t>n)  +  J  Bn{E' ,E)  |  E')dE' 


where 


H  |  *„(£))  =  E  |  9n(E)) 


It  satisfies  the  normalization  condition 

(*n(E)\*n(E'))=6(E-E') 

After  replacing  Eq.  58  into  Eqs.  59  and  60,  the  result  is 

I  Vn(En)  |2 


An{E)  r  = 


{E  —  En  —  An)2  +  7t2|  Vn(En)  |4 


=  p/ 


v.(g)  f 

E-E' 


dE‘ 


(58) 

(59) 

(60) 

(61) 

(62) 


Let  us  adopt  the  time-dependent  description  and  prepare  the  system  in  state  |  <j>n)  at 
t  =  0.  We  now  let  the  system  evolve  with  time  and  ask  what  is  the  probability  Pn(t)  of 
finding  the  system  still  in  state  |  4>n)  state  at  time  t?  The  answer  is 

Pn(t)  =  exp(—t/rn)  (63) 


where 

-  ft 
T"  ~  2ir\  Vn(En )  |2 

The  above  discussion  can  be  easily  generalized  into  the  case  in  which  one  bound  state 
is  coupled  with  many  unbound  states  j  E,m),  as  happens  when  the  final  predissociated 
system  can  be  characterized  by  the  set  of  quantum  numbers  m  describing  the  internal 
states  of  the  fragments.  The  result  is 


(64) 


_ ft _ 

**  Hm\V?(En)  I" 


(65). 


with  the  coupling  between  the  bound  state  |  <f>n)  and  the  mth  unbound  state  \  E,m)  as 


Vn{E)  =  {<t>n\H  \  E,m) 


(66) 
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5-2.2  The  Form  of  Coupling  Element  in  Nuclear  Coordinates 


We  now  discuss  the  nature  of  the  coupling  elements  which  appear  in  Fano's  formulas 
(Eqs.  57,  66).  The  Hamiltonian  of  a  molecular  system  has  the  form: 

H  =  He  +Tn  (67) 

where  electronic  He  is  the  Hamiltonian  (with  the  nuclear  variables  as  parameters)  which 
contains  the  Coulombic  interactions  between  all  charged  particles  (including  nuclei)  in 
the  molecule,  and  Tjv  is  the  nuclear  kinetic  energy  operator.  We  assume  that  the  spin- 
containing  terms  are  negligible  compared  with  the  Coulombic  terms  as  is  the  case  for  H^. 
We  can  define  a  set  of  electronic  basis  functions  <f>n( r,  R)  which  satisfy 

tfe^(r;R)  =  F(R)^(r;R)  (68) 

(4>n  I  Mr  =  W  (69) 

where  the  r  denotes  the  set  of  all  electronic  coordinates  and  R  the  set  of  all  nuclear 
coordinates,  n  is  the  set  of  quantum  numbers  which  describe  the  eigenstates  of  He.  In 
order  to  solve  the  total  Schrodinger  equation 

=£*(  r,R)  (70) 

We  expand  the  wavefunction  in  the  <f>n  basis  set: 

*  =  X>»(R)<Mr;R)  (7i) 


From  Eqs.  67-71,  we  get 

£  7Mxn(R)<Mr;R)}  =  £(£-F„(R))xn(R)<Mr,R)  (72) 

n  n 


In  Eq.  72,  the  nuclear  kinetic  operator  acts  upon  both  nuclear  Xn(R)  and  electronic 
<£n(r;R)  wavefunctions.  No  choice  of  reference  frame  and  coordinate  system  used  to  de¬ 
scribe  the  molecular  system  is  so  far  implied. 

In  the  laboratory  reference  frame,  the  nuclear  kinetic  energy  operator  Tn  is  given  by 


where  the  sum  is  over  all  nuclei  with  /z,  as  their  masses.  When  Eq.  73  is  put  into  Eq.  72 
with  the  help  of  Eqs.  68  and  69,  we  get  the  equations  which  Xn  satisfy: 

h2 

TsXn  +  ^ - 4>n^ i4>n'  '  ^ iXn'  +  ^  Xn'  =  {E  —  Vn)xn  (74) 

n',i  ***  n' 

In  the  normal  treatment  of  the  nuclear  motion,  only  one  term  is  retained  and  the 
coupled  nuclear  equations  are  replaced  by  a  set  of  uncoupled  equations 

r„xS°(R)  =  (E-V„(R))x®0(H.)  (75) 

The  corresponding  electro-nuclear  wavefunction  is  then 

¥°°(r,R)  =  xS°(R)A.(r;R)  (76) 


It  is  easy  to  show  that  the  matrix  element  of  the  total  Hamiltonian  H  between  a  pair 
of  wavefunctions  in  the  form  of  Eq.  76  but  with  different  electronic  states  is: 


-ft3 


!  H  I  =  (x®°  1 I  — v,  I  +  (*.  I  T„  I  *„.)  I 


x°'°> 


(77) 


In  our  application  to  predissociation,  the  initial  nuclear  bound  state  and  the  final  nuclear 
unbound  state  are  Born-Oppenheimer  solutions  (on  different  electronic  potential  energy 
surfaces),  and  the  corresponding  predissociation  coupling  matrix  elements  needed  for  ap¬ 
plication  of  the  Fano's  theory  are,  in  Cartesian  coordinates,  those  shown  in  Eq.  77. 


In  practice,  the  reference  frame  for  the  electronic  motion  is  usually  chosen  to  be  the 
body-fixed  frame  of  the  nuclei  and  the  translational  motion  of  the  center  of  mass  of  the 
system  is  removed.  Because  of  the  large  mass  difference  between  the  nuclei  and  electrons, 
the  center  of  mass  of  the  molecule  is  very  close  to  the  center  of  mass  of  the  nuclei,  and 
their  difference  can  be  safely  neglected  for  the  present  purposes.  The  coupling  elements 
appearing  in  Eq.  77  can  be  expressed  in  terms  of  elements  of  the  type 


(78) 


and 

(<k»  I  -q~2  I  (79)' 

where  x  can  be  one  of  the  internal  nuclear  hyperspherical  variables  (p,  u;*,  7>)  i.c.  the 
hyperradius  and  the  two  internal  hyperangles.  The  elements  ( 4>n  \  ^  |  <t>n')  and  (<pn  j  \ 
4>n’)  must  now  be  calculated  with  the  help  of  the  corresponding  electronic  wavefunctions. 
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6.  Summary 


In  order  to  assess  its  potential  as  a  future  rocket  propellant,  the  lifetime  and  de¬ 
cay  mechanisms  of  the  Hz  molecule  in  its  metastable  2 pz  2A^  electronic  state  must  be 
understood.  We  initiated  a  theoretical  study  of  this  decay  lifetime,  in  parallel  with  the 
experimental  work  on  the  properties  of  this  species. 

We  have  identified  the  necessary  steps  for  the  theoretical  investigation.  We  initiated 
the  calculation  of  the  Hz  electronic  potential  energy  surfaces  for  the  lowest  four  electronic 
states  (2 px<y  2E'{  1,2),  2s  2A\  and  2pa  2A.2#)  ,  and  of  the  electric  transition  dipole  moments 
among  them.  These  calculations  will  be  completed  in  the  near  future. 

We  have  developed  a  very  general  and  powerful  hyperspherical  coordinate  propagation 
method  to  obtain  the  ro-vibrational  nuclear  bound  states,  especially  useful  for  Az  systems 
(like  Hz)  having  three  identical  nuclei,  for  which  nuclear  permutation  symmetry  plays  an 
important  role.  We  also  used  for  such  calculations  the  variational  method  developed  by 
Tennyson  and  Sutcliffe  as  a  general  and  robust  way  to  treat  the  triatomic  ro-vibrational 
motion.  Both  methods  have  been  tested  on  the  upper  manifold  of  the  DMBE  surfaces22 
(2p  2E‘(2)),  with  or  without  inclusion  of  the  molecular  Aharanov-Bohm  effect68.  We  have 
performed  scattering  calculations  on  the  lower  manifold,  up  to  energies  of  2.4  eV,  using 
hyperspherical  coordinates62-64,69,  and  also  showed  how  to  use  parallel  computers  for  such 
calculations80. 

We  have  developed  a  numerical  methods  for  generating  scattering  wavefunctions  via 
backward  propagation,  once  the  scattering  matrix  has  been  calculated.  Such  wavefunctions 
are  needed  for  the  predissociation  lifetime  calculations  of  interest. 

After  the  ro-vibrational  wavefunctions  on  both  the  2s  2  A\  and  2 pa  2 A'£  potential 
surfaces,  and  the  electric  transition  dipole  moment  between  them  become  known,  the  ra¬ 
diative  lifetime  calculation  of  the  2 pz  2  A'2  state  to  the  2s  2  A\  state  can  be  finally  performed 
accurately.  Calculation  of  the  predissociative  lifetime  of  the  2 pa  2A'2  state  will  require  a 
calculation  of  the  appropriate  coupling  elements  to  the  ground  and  first  excited  electronic 
states. 
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An  efficient  numerical  method  of  calculating  surface  functions  for  accurate  quantum  mechanical  three-dimensional  reactive 
scattering  using  symmetrized  hyper-spherical  coordinates  has  been  developed.  This  method  is  at  least  20  times  faster  than  the 
finite-element  method  used  previously  and  us  accuracy  is  demonstrated  for  the  HrH.  system. 


1.  Introduction 

Accurate  quantum  solutions  for  three-dimensional  reactive  scattering  for  triatomic  systems  were  first  cal¬ 
culated  in  the  mid  1970s  for  the  system  H  +  H:  [1-5],  The  difficulty  and  computational  expense  of  these  cal¬ 
culations  has,  until  recently,  precluded  extension  to  higher  energies  and  more  complex  systems;  however,  the 
development  of  more  efficient  algorithms  coupled  with  increased  access  to  supercomputers  has  resulted  in  a 
resurgence  of  activity  in  this  field  [6-17],  In  particular,  the  use  of  symmetrized  hvper-spherical  coordinates 
(SHC)  and  local  hyper-sphencal  surface  functions  (LHSF)  [18,19]  is  a  very  promising  approach. 

The  first  accurate  calculations  of  3D  reactive  scattering  matrices  using  a  hyper-spherical  coordinate  method 
were  recently  performed  on  the  total  angular  momentum  J=Q  partial  wave  of  the  H  +  H2  system  [6],  This 
method,  applied  to  the  PK.2  potential  energy  surface  [20],  involved  the  calculation  of  sets  of  LHSF  using  a 
two-dimensional  finite-element  (FE)  approach.  The  FE  method  is  accurate  and  reliable  for  this  system,  and 
has  been  used  to  extend  the  range  of  energies  at  which  the  corresponding  7=0  partial  wave  scattering  matrices 
have  been  calculated  to  1.6  eV  [11];  however,  extension  to  higher  values  of  J  and  to  less  symmetric  systems 
requires  an  excessive  increase  in  computational  effort.  As  a  consequence,  there  was  a  need  to  develop  a  more 
efficient  method  for  calculating  these  LHSF. 

In  this  paper  we  present  a  new  variational  method  for  calculating  LHSF.  The  formalism  is  described  in  sec¬ 
tion  2.  Section  3  discusses  the  numerical  parameters  used  and  section  4  compares  the  results  of  LHSF  and 
scattering  calculations  for  the  7= 0  partial  wave  on  both  the  PK.2  and  LSTH  [21,22]  potential  energy  surfaces 
with  those  of  previous  calculations  using  a  finite-element  method  [6,1 1.23]  to  obtain  the  surface  functions. 
In  addition,  some  comparison  of  the  7=  1  PK2  scattering  results  with  those  of  the  matching  method  [3]  are 
made.  A  summary  is  given  in  section  5. 
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2.  Formalism  of  variational  surface  functions 


The  SHC  coordinate  system  used  in  this  calculation  has  been  described  previously  [6.1 1.18.19]  Let  A„.  A,,, 
A  be  the  atoms  of  a  tnatomtc  system  and  A,  v .  k  an  arbitrary  permutation  of  a,  /?.  y.  The  A  SHC  for  this  system 

are 


p=  (r\  +Rj)u2,  cJi  =  2  arctan  ,  y^sarccos 

A; 


Rrr. 


’A  $  n  , 


(1) 


where  rx  is  the  mass-scaled  [24,25]  intemuclear  vector  for  the  diatom  A, A,  and  Rx  the  position  vector  of  the 
atom  Ax  with  respect  to  the  center  of  mass  of  AA*-  The  orientation  of  the  system  in  space  is  determined  by 
the  Euler  angles  8;,  <?x  (the  polar  angles  of  Rx  with  respect  to  a  space-fixed  OZ  axis ),  and  ip-  (the  angle  between 
the  Rx,rx  and  Rx,OZ  half-planes).  In  this  coordinate  system  the  Hamiltonian  is  expressed  as 


,5/2. 


A- 
2  UP2 


+  V(p,  OJX,  '/;  )  + 


\5fr_ 

8 UP2  ’ 


(2) 


in  which  the  global  reduced  mass  p  is  defined  as  \mxmumKl (mx  +  w„+  mK)  ]|,:.  The  generalized  or  grand  ca¬ 
nonical  angular  momentum  operator  A 2  is  defined  by 

A2  =  L;+ - - -  +  .  Jfr - -,  (3) 

cos  *  ( 1  oj, )  sin-  ( ( wj 

where  I l  and  f~  are  the  angular  momentum  operators  associated  with  the  vectors  R.  and  r „.  respectively,  and 


Li 


—  4fi~  (  6- 
sin  <dx  Veto; 


sin  (jjx 


(4) 


is  an  angular  momentum  associated  with  the  hyper-angle  w;.  The  term  V(p.  ojx,  yx)  is  the  potential  energy 
function  of  the  electronically  adiabatic  tnatomic  system. 

The  equation  that  defines  the  LHSF  <pJ„K,nr  with  associated  eigenvalues  tinr  is 


+  V(p,  ojx,  y,  )'j<t>J"nr(L\p)=tinr(p)  <PiMnr(^:p) , 


(5) 


in  which  Ci  stands  for  the  set  of  five  hyper-spherical  angles  (ai*,  yx,  9X,  <t>x,  Vi)-  The  indices  J,  M,  17,  T  are, 
respectively,  the  quantum  numbers  of  the  total  angular  momentum  of  the  system,  its  projection  on  the  space- 
fixed  OZ  axis,  the  inversion  parity  of  the  triatomic  system  through  its  center  of  mass,  and  the  irreducible  rep¬ 
resentation  of  the  surface  function  in  the  permutation  group  of  the  system  (P3  for  H  +  H:).  The  boundary  con¬ 
ditions  for  eq.  (5)  are  the  usual  “well-behavedness'  ones  (single  valuedndss.  continuity,  non-divergence, 
differentiability,  etc. ).  The  index  n  denotes  a  quantum  number  which,  in  addition  to  J.  \f,  17.  F.  uniquely  labels 
the  LHSF.  A  set  of  two-dimensional  surface  functions  viificv,,  >'*;/?)  independent  of  the  orientation  of  the 
s”  ’em  in  space  can  be  defined  by  expansion  of  the  LHSF  in  terms  of  Wigner  rotation  matrices 
dx,  V, )  [26]: 


r(L\P)=  z  ZiTaid), ,  e, ,  V, )  H/inar(  <J..,yA\p) 
a- 0 


(6) 


where 

9J^o(<t>x,  V,)  =  A'jo[Dito(0,,d;,  <px)  +  (-\  )J*n*aDif..o(0i,  6i,  Vi)\ 

and  Sja  is  a  normalization  constant.  is  even  (odd)  with  respect  to  inversion  of  the  system  through  its 
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center  of  mass  for  77=0  ( 1 ),  and  in  the  case  of  72  =  0  is  zero  when  J+Tl  is  odd.  The  boundary  conditions  for 
the  if/iar ( co,,  p)  which  result  from  the  “well-behavedness"  of  the  (P'n',nr  are  that 


vri2r(c«Ji.>’i  =  {0.  rt ):p)  =0 

and 


for  72*0 


t —  [r!£>r(Mi'7i:p)]\n-io  .;=0  for 72=0  . 
oy< 

Furthermore,  the  potential  function  V(p,  cox,  yx)  has  in  general  extrema  at  y„  =  {0,  it}  (corresponding  to  col- 
linear  configurations  of  the  system).  As  a  result  of  these  considerations,  the  yiar  can  be  factored  as  a  function 
of  cUi  times  a  normalized  (2)  associated  Legendre  function  ^f(cosyx)  in  the  vicinity  of  y^  =  0  and  n.  This 
makes  it  both  convenient  and  desirable  (because  of  the  presence  in  eq.  ( 3 )  of  the  operator  jx )  to  expand  these 
according  to 

Winar(^.7,-.p)=  £  (cosy„)  <~~~  ■  <7> 

,zh  sinoi; 

where  the  coefficients  \p)  are  called  one-dimensional  surface  functions.  Replacement  of  eqs.  (  7 )  and 

(6)  into  eq.  (5)  leads  to  the  equation  satisfied  by  these  functions: 

ff  J  .  ,v  j(j+\)+ju+\)-2a1 .  yo+i)  1  - - 


HHJi 


COS:(  Jcu;  ) 


70+ 1)  1 
tn:(  lcj;)  J 


YT—  O)  {.O'.  Q)  0i%,t((oi-,p)+i.{J.  Q)  {-O'.  G)  0i%-uicot:p)  ] 


_ fr_ 

COS2(  id);  ) 

f  v  afip.  OJ; )  oi% ( <x)A:  p )  =  t inr(p )  0J„% (aip.p) 

f-a 


The  multiplicative  factor  ( sin  cu,  ) _ '  has  been  introduced  into  eq.  (7)  because  of  the  form  of  eq.  (4).  The 
presence  of  this  term  forces  the  boundary  conditions  <p^(w;  =  {0.  n};  p)  =0  for  eq.  (8).  These  conditions  are 
necessary  :or  not  to  diverge  at  01,=  [0,  ji)  but  may  not  be  sufficient;  however,  in  practice  they  have  indeed 
sufficed  for  H  +  H;,  and  we  do  not  anticipate  problems  with  other  systems.  In  eq.  (8),  {.(7.  k)  = 
[7(7+  1  )-k(kz  1 ) ] l/:  and  the  term  V%  is  given  by 
* 

V°Ap  ,Wi  -  I  (  cos  y, )  V(p.  10, Jf(  cos  yx )  sin  >•„  dy>  .  (9) 


It  is  impu  nant  to  note  that  for  eq.  ( 8  )  to  be  valid,  the  functions  with  72=0  must  be  defined  to  be  iden¬ 
tically  equal  to  zero  when  J+ 77  is  odd.  The  set  of  equations  ( 6  ).  (7).  and  ( 8 )  are  equivalent  to  eq.  ( 5  ). 

The  vr  lationai  basis  set  is  suggested  by  the  expansion  equation  (7)  and  by  eq.  (8).  We  define  functions 
tff0((oy,  f  I  with  associated  eigenvalues  eJ,,a(p )  which  satisfy  the  latter  after  the  72  and  j  coupling  is  removed: 


'  _  1*11  i_  + 

j(j+\)+j(j+\)-2a2 

70+1) 

pp:\  ai: 

4cos;(  (coj 

4sin*(Jcu;)/  "  J 

These  fui  :tions  are  required  to  satisfy  the  same  boundary  conditions  as  the  iua(0\p)  =  t^n(KP)  =  0. 

Wc  now  oefine  a  five-dimensional  vari.uional  baso  set  by 

FJ”!]{Z;.p)  =  7/{7n(OA.0..v, )  .j»f'(cosyJ/^(W;.p)  .  (11) 

where  for  notational  convenience  we  have  defined /^(oj, ; p)  =  w;: p ) /sin  Since  this  basis  set  is  con- 
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centrated  in  the  /.  arrangement  channel  region  when  p  is  sufficiently  large,  accurate  representation  of  the  surface 
functions  concentrated  in  the  v  and  k  channels  may  require  a  large  number  of  terms  at  such  values  of  p.  To 
overcome  this  difficulty,  a  new  basis  set  is  constructed  which  consists  of  the  union  of  the  basis  sets  F{''JJ . 

and  Fi'lfi .  Furthermore,  for  systems  containing  either  two  or  three  identical  atoms  we  construct  sym¬ 
metrized  basis  sets  which  belong  to  irreducible  representations  Tof  the  P;  or  Pj  permutation  groups,  respec¬ 
tively.  For  the  three  identical  atom  case,  these  symmetrized  variational  basis  sets  are  given  by 

FJ.$nr(H-,p)=  1  cfy.Ff'2<C, ;/»  .  (12) 

A 

where  the  sum  in  X’  is  over  X,  v  and  k,  the  cf,  are  easily  determined  constants,  and  the  sets  of  angles  £,  and 
C,  are  considered  to  be  functions  of  Cj-  The  functions  F/ya/7r  will  be  referred  to  as  primitives  to  distinguish 
them  from  the  unsymmetrized  basis  functions  F{"3  .  The  five-dimensional  LHSF  are  now  expanded  in  terms 
of  these  primitives: 

<*>  i',nr(  &:/»«!  aUS. (P )  nn  C -P >  ■  (13) 

•/a 

The  primitive  basis  set  is  not  orthogonal,  since  the  variational  basis  sets  with  different  overlap:  therefore, 
calculation  of  the  coefficients  requires  the  determination  of  overlap  integrals  for  the  variational  basis  set 
as  well  as  integrals  involving  the  Hamiltonian.  Integration  over  the  three  Euler  angles  6,,  0,.  and  <p,  is  analy  tic, 
leaving  two-dimensional  quadratures  to  be  done  numerically.  These  quadratures  are  the  most  expensive  part 
of  the  entire  computation.  Any  quadrature  scheme  may  be  employed:  the  one  we  used  is  discussed  in  section 
3. 

Once  all  of  the  necessary  integrals  have  been  calculated,  the  a(%f„  coefficients  are  determined  by  a  generalized 
eigenvalue-eigenvector  procedure.  With  sufficienty  large  basis  sets,  the  overlap  matrix  between  the  primitives 
becomes  nearly  singular  as  a  consequence  of  near  linear  dependence:  for  this  reason  it  was  necessary  to  develop 
a  method  for  dealing  with  this  situation.  The  primitive  overlap  matrix  is  diagonalized,  and  eigenvectors  cor¬ 
responding  to  eigenvalues  smaller  than  a  tolerance  parameter  (for  the  calculations  described  below,  this  pa¬ 
rameter  was  set  to  zero )  are  eliminated  from  consideration  The  remaining  set  of  eigenvectors  is  used  to  transform 
the  Hamiltonian  matrix  to  yield  a  new  eigenvalue  problem  from  which  the  linear  dependence  effects  have  been 
removed. 

From  eqs.  (6),  (12),  and  (13),  the  expansion  of  the  two-dimensional  surface  functions  vJ„ar  yp  p)  in 
terms  of  the  one-dimensional  functions /"a  can  be  shown  to  be 

VJnar(w>.'Vi\P)=  I  a{?a*{c&Q  [•  +  (-!  )J~n6$)??  (cos*a)/"0  (0J,\p) 

„a 

+  cr„dJs?a  (d„J  9°  (cos y,)fija  (cu„;p)+c£,(  - 1  )n*°d^  (d„„)  (cos;v)/»o  (w„;/>)  .  (14) 

The  functions  daa  (J)  —dJan  (d)  +  (  —  1  )J*n*a  dJa  (d).  where  dJ  is  the  Wigner  little  d  matrix  (2-»):  they 
appear  in  eq.  (14)  because  of  the  integration  of  products  of  two  functions  depending  on  different  Euler 
angles.  The  angles  J„,  between  the  vectors  R,  and  Rv  and  d,„  between  the  vectors  RA  and  Rx  are  functions  of 
oji  and  yA  only. 

The  same  formalism  is  used  in  determining  the  surface  functions  at  values  of  p  for  which  the  surface  function 
amplitude  is  negligible  in  regions  of  configuration  space  in  the  interstices  between  the  arrangement  channel 
regions.  For  such  values  of  p,  the  overlap  between /  Jv‘a  and/ /'  o  vanishes,  making  the  set  of  primitive  func¬ 
tions  be  automatically  orthogonal;  this  greatly  reduces  the  numerical  work  necessary  for  the  LHSF  calculation, 
because  the  basis  set  includes  only  the  X  basis  functions. 

The  calculations  of  the  six-dimensional  scattering  wavefunction  is  done  by  expanding  it  in  terms  of 

the  five-dimensional  LHSF: 
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(15) 


The  <Pi',,,r  are  determined  at  a  set  of  discrete  values  of  t\  labeled  pr  Substituting  eq.  (15)  into  the  (time- 
independent)  Schrodinger  equation  corresponding  to  the  Hamiltonian  defined  by  eq.  (2)  and  using  eq.  (5). 
the  coefficients  b(p\p,)  are  found  to  satisfy 


V  2ppir‘ 


d: 


d P 


iP 


15ft* 
8  HP2 


■©' 


tJnnr(p,  )-E)K  binr(p\  p,)+  X  binr(p\  p,)[J  jnr] "  (p\p,)=0  ,•  (16) 


in  which  the  interaction  matrix  Jjnr  is  defined  by 
l JjnrVn  (p. P,)  =  <<p i',nr( C ; p, ) I  f(p, OJ, , 7i \p,)\<?> n,nr( C ; p, ) > 

=  X  I  ai}arn(p,)ai%  AP^<Ft}anr (Ci\P,)\i',(P- toi,Yi,P,)\F J,yzr {&.(>,) >  ,  (17) 

i  jU  v  j  Q 

with  P(p,  oi;.,  y,  -P,)-  l’(p.  cu;  ,  •/; )  -  (p,/p ) 2  V(p, ,  Cii^,  ft ).  The  integrals  in  the  right-most  part  of  eq.  ( 17 )  are 
obtained  from  linear  combinations  of  related  integrals  involving  the  variational  basis  set  (11). 

The  coefficients  b(p:p;)  are  calculated  as  a  function  of p  in  a  region  near p,  corresponding  to  a  hyper-spher¬ 
ical  shell.  The  smooth  matching  of  the  scattering  wavefunc'  on  across  the  boundary  ,  of  adjacent  hyper- 
spherical  shells  is  accomplished  by  imposing  the  condition" 


^(p,.,.  ,;/?'+,)  =  X  binr(p,,  +  ,:p,)[(Ljnr):  (a-i.A)  . 


(18) 


=  ^  f?bfjp-.p,)\  [CJnr]:(-+i  -i} 
\  *P  A-p'.V,  -  V  bp 


in  which  the  overlap  matrices  Cmr  are  defined  by 
[Cjnr]l  (p’.i.p~)  =  (<t>J„'mr(AP,~\)\<& iunr( Ct : P, ) > 

=  X  I  ailhrAp,,i)a{%Ap,)(Fi%nr(Zi\pl+>)\Fi”2r(Zi-,p,)'>. 


i  /Q  i  j  Q 


(19) 


(20) 


The  methodology  described  above  is  closely  related  in  spirit  to  the  method  independently  developed  by  Schatz 
[17],  which  was  published  after  the  present  work  was  completed.  The  major  differences  are  in  the  selection 
of  reference  potential  for  calculation  of  the  c/ft-dependent  portion  of  the  basis  set  and  in  the  method  for  dealing 
with  overcompleteness  of  the  basis  set.  Our  reference  potential,  denoted  by  the  term  V„  in  eq.  (10),  is  the 
potential  energy  surface  at  fixed  p  averaged  over  the  diatomic  rotation:  the  choice  of  this  reference  potential 
naturally  follows  from  the  expression  for  the  one-dimensional  surface  functions  0J„%f  of  eq.  ( 8 ).  We  allow  for 
the  large  amount  of  linear  dependence  which  is  produced  by  this  reference  potential  at  small  values  of  p  by 
the  method  for  solution  of  the  generalized  eigenvalue  problem  described  above.  This  does  not  increase  the  time 
required  for  the  calculation.  Schatz.  on  the  other  hand,  chooses  for  his  reference  potential  V(p,  <y; ,  ft  = 
it./ 2 )  for  3.3  bohr  and  Hp  =  3.3  bohr.  to,_,  ft  =  it/ 2 )  for  p<  3.3  bohr;  the  change  in  the  reference  potential 
at  small  p  avoids  problems  with  linear  dependence. 


3.  Numerical  parameters 

One  of  the  most  important  parameters  in  the  calculations  performed  is  the  number  of  primitives  used  to 
expand  the  surface  functions.  In  addition  to  the  indices/.  M.  /7and  /"which  label  the  LHSF.  the  basis  functions 
and  the  primitives  formed  from  them  are  labelled  by  indices  v.j,  and  32.  which  asymptotically  correspond  re- 
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spectively  to  the  diatom  vibrational,  total  rotational  and  helicity  rotational  projection  quantum  numbers  [2.3  ]. 
The  range  of  (r ,j.  Q)  included  in  the  calculation  of  a  desired  set  of  LHSF  is  selected  by  preliminary  calculations 
at  a  small  subset  of  the  values  of  p  to  be  used  in  the  full  calculation.  For  these  initial  calculations,  the  basis 
set  is  deliberately  chosen  to  be  larger  than  necessary  for  the  accuracy  and  number  of  surface  functions  desired. 
The  number  of  accurate  surface  functions  obtained  by  this  method  is  determined  by  comparison  of  two  such 
calculations  with  different  size  basis  sets  at  eacft  value  of  p.  Examination  of  the  coefficients  of  the  basis  func¬ 
tions  contributing  to  each  surface  function-  considered  allows  selection  of  a  smaller  .basis  set  which  can  be  used 
with  minimal  loss  of  accuracy.  This  method  becomes  complicated  with  larger  variational  basis  sets,  due  to  in¬ 
creasing  linear  dependence  among  some  of  the  primitive  functions;  however,  the  overall  pattern  of  important 
coefficients  is  still  effective  in  optimizing  the  choice  of  basis  sets  for  succeeding  calculations.  When  this  method 
of  selection  is  used,  the  number  of  good  surface  functions  of  each  symmetry  which  are  produced  is  approxi¬ 
mately  one  half  the  number  of  primitives  of  that  symmetry  used  in  the  calculation. 

To  obtain  the  results  presented  below,  the  one-dimensional  numerical  functions  i^Q(uix\p)  from  eq.  (10) 
are  calculated  on  a  grid  of  450  ojx  points  using  a  one-dimensional  finite  element  method.  Each  element  is  qua¬ 
dratic  and  uses  two  Gauss-Legendre  points.  The  reference  potential  V°  for  these  functions  is  determined  by 
a  Gauss  -Legendre  quadrature  with  96  points.  The  grid  for  the  two-dimensional  integrals  is  the  direct  product 
of  these  two  independent  quadrature  grids.  Convergence  of  the  surface  function  energies  with  the  fineness  of 
the  mesh  is  to  four  decimal  places. 

For  the  system  H  +  H;,  the  dependence  of  the  basis  functions  eq.  (11)  in  the  v  and  k  coordinates  is  the  same 
as  that  for  the X  channel  functions,  so  it  is  not  necessary  to  repeat  the  calculation  for  tifo  and  rl'a.  In  addition, 
the  integrals  between  products  of  functions  in  the  a  and  v  channels  equal  the  integrals  between  the  corre¬ 
sponding  functions  in  the  v  and  k  channels,  so  the  integration  need  only  be  done  for  X,v  pairs  to  obtain  the 
overlap  integrals  for  all  three  regions. 

LHSF  are  calculated  every  0.2  bohr  from  p  =  2.0  to  12.0  bohr.  and  interaction  matrices  Jjnr  (see  eqs.  (16) 
and  ( 1 7 ) )  are  determined  at  five  evenly  spaced  values  of  p  for  every  value  of  p,.  One  overlap  matrix  Cjnr  ( see 
eqs.  ( 18)- (20) )  is  calculated  between  sets  of  the  LHSF  at  each  pair  of  adjacent  p  values.  For  6.2  bohr  the 
variational  basis  functions  in  arrangement  channel  X  are  orthogonal  to  those  in  v  and  k  because  we  set  the 
maximum  value  of  cu;  equal  to  the  physically  reasonable  value  2  arcsin( 3.0 bohr/ p),  which  at  p=  6.2  bohr  equals 
57.9°  This  value  of  (a„  is  deep  in  a  classically  forbidden  region  for  all  y*  for  the  total  energies  discussed  below. 
As  a  result,  the  time  needed  for  the  surface  function  calculation  in  this  region  is  small  compared  to  that  for 
the  p<  6.2  bohr  region.  The  initial  value  problem  described  by  eq.  (16)  is  solved  using  a  logarithmic  derivative 
propagator  [27]  with  a  step  size  of  Ap=0.025  bohr  and  a  constant-p  projection  [6,1 1,28]  at  12.2  bohr.  These 
parameters  were  chosen  to  achieve  a  calculation  accuracy  about  equal  to  that  described  previously  [6.1 1 }. 

The  five-dimensional  basis  functions  are  generated  from  the  functions / according  to  eq.  ( 1 1 ).  For  each 
of  the  potential  energy  surfaces  and  for  7=0,  a  set  of  (v,  j,  J2=0)  quantum  numbers  was  chosen  to  give  a 
variational  basis  set  of  152  functions.  This  set  has  a  maximum  of  12  vibrational  functions  for  the  value  j  —  0. 
with  monotonicallv  decreasing  number  of  vibrations  for  each  succeeding  value  of  j  to  the  maximum  of  y  =  23 
for  which  only  one  vibrational  function  is  used.  Symmetrization  of  this  basis  set  yielded  76  A,.  76  A:  and  1 52 
E  primitives. 

For  the  LSTH  potential  energy  surface,  the  scattering  results  were  obtained  from  36  A,.  35  A;,  and  69  E 
LHSF  at  each  value  of  p.  The  calculation  of  each  LHSF  (including  the  evaluation  of  all  the  associated  overlap 
and  interaction  matrices  for  the  solution  of  the  propagation  equation  ( 16) )  requires  an  average  of  0.27  s  on 
a  Cray  X-MP/48.  The  timings  are  very  similar  for  the  PK2  potential  energy  surface. 

For  7=0,  the  variational  LHSF  calculation  is  about  a  factor  of  20  faster  than  the  finite-element  method  one 
for  equivalent  accuracy  [6,1 1  ].  We  estimate  that  the  numerical  effort  required  for  the  finite-element  calcu¬ 
lation  of  the  LHSF  will  increase  with  7  as  (7+  1 )“  with  2<a$3.  whereas  for  the  variational  method  a=s2; 
therefore,  the  speed  of  the  variational  LHSF  calculation  with  respect  to  that  of  the  finite-element  one  is  ex¬ 
pected  to  increase  with  increasing  7. 
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Fig.  1.  Probabilities  (a)  and  probability  differences  (b)  as  a 
function  of  total  energy  £  (lower  abscissa)  and  initial  relative 
translational  energy  £oo  (upper  abscissa)  for  the  /» 0  (0,  0, 
0 )  —  (0.  0.  0 )  E  symmetry  transition  in  H  +  H2  collisions  on  the 
PK2  potential  energy  surface.  The  symbol  (v,  j,  Q)  labels  an 
asymptotic  state  of  the  H  +  H2  system  in  which  v,  j  and  Q  are  the 
quantum  numbers  of  the  initial  or  final  H2  states  as  defined  in 
the  text.  The  vertical  arrows  on  the  upper  abscissa  denote  the 
energies  at  which  the  corresponding  H2(v./)  states  open  up.  The 
length  of  those  arrows  decreases  as  v  spans  the  values  0.  I,  and  2. 
and  the  numbers  0,  5.  and  1 0  associated  wuh  the  arrows  define  a 
labelling  for  the  value  of  j.  The  square  symbols  in  (a)  are  the 
current  variational  surface  function  results  and  the  solid  line  are 
the  FE  results  (11).  The  differences  between  the  former  and  the 
latter  are  plotted  in  (b). 


Eoc/eV 

00  02  04  06  08  10  12 


Fig.  2.  Same  as  for  fig.  I  for  the  LSTH  potential  energy  surface. 
The  FE  results  are  taken  from  ref.  [22 ). 


4.  Results  and  discussion 

As  can  be  seen  in  table  1,  the  present  variational  (V)  LHSF  energies  consistently  fall  below  those  calculated 
by  the  finite-element  method  ( FE )  [11,23],  with  a  maximum  reduction  of  about  65  meV  for  the  higher  energy 
LHSF.  As  both  methods  obey  a  minimum  principle,  this  implies  better  quality  of  the  LHSF  in  the  current 
method. 
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Flux  is  conserved  for  the  PK.2  surface  to  better  than  1%  for  energies  less  than  1.55  eV  for  the  scattering  cal¬ 
culations  using  FE  LHSF  [11]  and  for  energies  below  1.74  eV  for  the  current  calculations.  For  the  LSTH  sur¬ 
face.  flux  is  conserved  to  better  than  2°o  below  1.55  eV  for  the  FE  results  [23]  and.  for  the  current  results,  to 
better  than  l°o  below  1.68  eV  and  to  better  than  2°'o  between  1.68  and  1.74  eV.  Examination  of  the  scattering 
matrices  produced  by  each  method  shows  good  agreement  between  the  two  below  the  first  resonance  at  0.97 
eV,  with  a  difference  which  is  usually  no  greater  than  2%  for  probabilities  greater  than  10_:.  Above  this  energy 
the  agreement  between  the  results  for  the  v=0  to  v  =0  state  transitions  remains  equally  good;  however,  the 
agreement  for  v=  1  to  v'  =  1  transitions  is  not  as  good,  with  a  difference  usually  no  greater  than  4%  in  the  prob¬ 
abilities  greater  than  10_:  for  these  transitions.  On  the  basis  of  the  lower  LHSF  energies  and  the  better  scat¬ 
tering  matrix  unitarities  of  the  current  method,  we  believe  that  the  current  scattering  calculations  are  more 
accurate  than  the  ones  using  FE  LHSF.  A  comparison  of  7=  0  probability  curves  generated  by  the  two  methods 
is  plotted  in  fig.  1  for  the  PK2  surface  and  in  fig.  2  for  the  LSTH  surface,  together  with  the  difference  of  the 
results  of  the  two  methods.  The  relative  differences  for  the  probabilities  of  these  figures  never  exceed  6%  for 
the  PK2  or  2%  for  the  LSTH  potential  energy  surface;  the  greater  maximum  relative  difference  for  the  former 
is  due  to  the  smaller  minimum  value  of  the  probability  itself  rather  than  to  a  greater  difference  between  the 
results. 

The  scattering  results  calculated  from  the  variational  LHSF  are  converged  with  respect  to  number  of  surface 
functions  used  in  the  propagation  to  2%  in  the  probabilities  greater  than  0.1  and  1.5 5  in  the  corresponding 
scattering  matrix  element  phases.  The  energies  of  the  LHSF  corresponding  to  asymptotically  open  states  are 
converged  to  0.5°o  with  respect  to  size  of  the  basis  set  used  in  their  calculation,  and  thus  the  scattering  results 
are  also  well  converged  with  respect  to  this  parameter. 

Calculations  were  also  performed  for  both  parities  of  the  7=  1  partial  wave  and  are  of  similar  quality.  7=  1 
probabilities  in  excess  of  0.1  generally  agree  with  the  matching  method  ones  [3]  (which  were  obtained  up  to 
0.7  eV  only)  to  better  than  about  3%.  In  none  of  these  calculations  have  we  encountered  the  difficulties  pre¬ 
viously  predicted  [10].  A  detailed  analysis  of  the  7=  1  results,  up  to  1.75  eV.  and  of  the  corresponding  res¬ 
onances  will  be  the  subject  of  a  separate  publication. 


5.  Summary 

A  new  general  variational  method  for  calculating  local  LHSF  was  described.  It  is  about  a  factor  of  20  more 
efficient  than  the  finite-element  method  for  the  7=0  partial  wave  of  H  +  H:;  this  relative  efficiency  is  expected 
to  increase  with  increasing  7.  The  results  of  the  two  methods  agree  well  at  the  LHSF  level,  and  at  the  scattering 
matrix  level  agree  well  for  energies  below  0.97  eV  and  moderately  well  for  higher  energies;  the  variational  ones 
are  believed  to  be  the  more  accurate  ones. 
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We  have  performed  accurate  quantum  mechanical  three-dimensional  reactive  scattering  calculations  for  both  panties  of  the 
J-  1  partial  wave  of  the  H  +  H;  system  up  to  total  energies  of  1.75  eV.  The  collision  lifetime  resonance  spectra  for  both  y=0  and 
J=  I  are  discussed  in  terms  of  the  characteristics  of  the  system's  potential  energy  surface  and  of  a  simple  physical  model  involving 
its  symmetry  properties. 


I.  Introduction 

We  have  recently  developed  a  new  variational 
technique  [  1  ]  for  calculating  the  local  hyperspher- 
ical  surface  functions  (LHSF)  necessary  for  per¬ 
forming  three-dimensional  (3D)  quantum  mechan¬ 
ical  reactive  scattering  calculations  by  the 
symmetrized  hyperspherical  coordinate  method.  We 
have  shown  that  this  technique  produces  results  of 
similar  quality  as  the  finite  element  (FE)  one  pre¬ 
viously  used  12.3],  with  significantly  less  numerical 
effort.  Using  the  LHSF  generated  by  this  variational 
method,  we  have  performed  3D  reactive  scattering 
calculations  for  the  H  +  H-  system  on  the  PK2  [4] 
and  LSTH  [5.6]  potential  energy  surfaces  for  the 
7=0  and  both  parities  of  the  7=  1  partial  waves.  The 
calculations  are  of  sufficiently  high  quality  for  res¬ 
onance  analysis  using  the  collision  lifetime  matrix 
formalism  [7],  We  briefly  describe  the  parameters 
used  in  the  calculation,  and  follow  with  a  presenta¬ 
tion  and  analysis  T  the  results.  Other  recent  hyper- 
sphertcal  ca'culations  for  J=\  H  +  H-  have  been 
published  by  Schatz  [  8  ]  and  by  Pack.  Parker  and  co- 

1  Work  performed  in  partial  fulfillment  of  the  requirements  for 
the  Ph  D.  Degree  in  Chemistry  at  the  California  Institute  of 
Technology 

;  Current  address.  Mail  code  206-49.  California  Institute  of 
Technology.  Pasadena.  CA  9 1  1 25.  USA. 
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workers  [9],  over  a  more  limited  energy  range.  In 
addition,  results  using  different  methods  have  been 
obtained  by  Mladenovic  et  al.  [10]  and  by  Zang  and 
Miller  [II], 


2.  Method 

The  details  of  the  partial  wave  methodology  for 
arbitrary  /  and  the  parameters  used  for  the  7=0  cal¬ 
culation  have  been  presented  previously  [  1  ] ;  there¬ 
fore,  only  a  few  relevant  points  will  be  mentioned 
here.  The  7=0  and  7=  1  calculations  use  the  same 
values  for  many  numerical  parameters:  the  choice  of 
gnd  is  the  same,  as  are  the  number  and  location  of 
the  sets  of  LHSF.  The  basis  set  for  7=0  is  formed 
from  a  choice  of  quantum  numbers  (r.  j,  Q= 0) 
[1.12],  yielding  a-total  of  152  functions.  This  basis 
is  symmetrized  according  to  the  irreducible  repre¬ 
sentations  of  the  P,  symmetry  group  of  H  +  H:.  to 
give  76  A,,  76  A;.  and  152  E  “primitive"  functions. 
These  are  used  as  a  variational  basis  for  calculation 
of  the  LHSF  of  the  corresponding  symmetry.  The 
same  set  of  r  and  j  quantum  numbers  is  used  for  the 
7=  1  calculations:  in  addition,  making  equal  to  both 
0  and  1  produces  a  variational  basis  set  of  292  func¬ 
tions.  From  these,  a  primitive  basis  set  is  generated 
consisting  of  64  A,,  76  A,  and  140  E  primitives  for 
the  even  parity  77=0  (which  for  J=  1  contains  only 
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72  =  1  functions)  and  another  with  140  A„  152  A: 
and  292  E  primitives  for  the  odd  parity  77=  1  (con¬ 
taining  both  72  =  0  and  72=  1  functions). 

The  7=  1  scattering  results  on  the  P1C2  surface  were 
obtained  from  31  A,.  31  A:,  and  64  E  LHSF  for  77=0 
and  67  A„  62  A;,  and  133  E  LHSF  for  77=  1.  The 
calculation  of  each  7=  1  LHSF  and  all  associated 
matrices  used  in  the  logarithmic  derivative  propa¬ 
gation  [13]  required  an  average  of  1 3. 1  s  on  an  SCS- 
40  minisupercomputer,  as  compared  to  6.9  s  for  7=0. 
Similarly,  the  7=  1  scattering  calculations  for  the 
LSTH  surface  used  32  A,,  32  A2  and  64  E  LHST  for 
77=0  and  74  A„  70  A:,  and  127  E  LHSF  for  77=  1, 
with  an  average  time  of  12.5  s  per  LHSF,  compared 
to  6.6  s  for  7=  0.  The  corresponding  maximum  de¬ 
viation  from  flux  conservation  is  never  greater  than 
1%  for  the  PK2  and  2%  for  the  LSTH  surface. 


3.  Results  and  discussion 

From  our  irreducible  representation  scattering 
matrices  we  have  calculated  distinguishable  atom 
7=  1  state-to-state  reaction  probabilities  over  the  en¬ 
ergy  range  0.3  to  1.73  eV  for  the  PK2  potential  en¬ 
ergy  surface.  The  coupled  channel  (CC)  results  pub¬ 
lished  by  Schatz  at  0.5  and  0.6  eV  [  8  ]  and  ours  agree 
to  within  10%,  which  is  reasonable  since  he  used  a 
much  smaiier  basis  set  than  ours.  Schatz  also  made 
calculations  based  on  the  coupled  states  (CS)  ap¬ 
proximation  using  a  larger  basis  set  than  in  this  CC 
method.  These  CS  probabilities  are  closer  to  our 
highly  converged  values  than  the  CC  ones,  indicat¬ 
ing  that  the  CS  approximation  for  his  larger  basis  set 
is  more  accurate  than  the  CC  results  using  his  smaller 
basis  set.  Our  7=  1  results  agree  with  the  LSTH  cal¬ 
culations  of  Zhang  and  Miller  [11]  at  1.14  eV  to 
about  2%  or  better. 

We  performed  lifetime  matrix  analysis  [  3.7. 1 4  ]  of 
each  of  the  matrices  Sjnr  for  T'=A,,  A2,  E,  77=0,  1, 
and  7=0, 1.  We  label  the  resonances  obtained  by  the 
notation  appropriate  for  vibrational  states  of  linear 
tnatomic  molecules.  ( vu  p*.  r3),  where  p,.  p2,  and  p3 
denote  respectively  the  quantum  numbers  for  the 
symmetric,  bend  and  asymmetric  vibrations  and  K 
is  the  quantum  number  of  the  vibrational  angular 


momentum  [15]  .  For  the  resonance  state,  r,.  r; 
and  ty  denote  approximate  constants  of  the  motion; 
their  values  are  chosen  on  the  basis  of  the  energy 
spacing  of  the  resonances.  In  the  present  paper,  as 
well  as  in  all  previous  ones  using  this  labeling,  it  has 
been  customary  to  set  p3  =  0  [3,16.17  ],  implying  that 
such  resonances  have  no  asymmetric  stretch  char¬ 
acter.  However,  modeling  of  collinear  H}  resonances 
(for  which  only  v,  and  v}  are  defined)  has  shown  that 
they  may  have  significant  asymmetric  as  .well  as 
symmetric  stretch  character..  The  vibrationally  adi¬ 
abatic  model  suggests  that  the  lowest  collinear  H3 
resonance  be  assigned  the  quantum  numbers  p,  =  1 . 
p3  =  0  [18],  whereas  the  hvperspherically  adiabatic 
model  leads  to  the  assignment  p,  =0,  p3  =  2  [  19,20], 
corresponding  to  the  second  excited  state  of  the 
asymmetric  stretch  and  asymptotically  correlating  to 
the  p=  1  state  of  the  isolated  diatom.  Therefore,  the 
nodal  structures  of  the  corresponding  model  wave- 
functions  are  completely  different,  and  neither  should 
be  assumed  correct  without  further  comparison  with 
the  accurate  resonance  wavefunction.  The  assign¬ 
ment  p3  =  0  used  in  this  paper  corresponds  to  a  vi¬ 
brationally  adiabatic  description,  but  is  a  matter  of 
notation  rather  than  of  physical  validity. 

Lifetime  matrix  analyses  of  the  7=0  scattering 
matrices  for  the  PK2  surface  were  previously  per¬ 
formed  up  tc  1.6  eV  using  the  FE  method  for  cal¬ 
culating  t  ie  LHSF  [  3  ] .  They  were  recalculated  using 
the  variational  LHSF  approach,  and  the  results  are 
comparable.  The  resonant  time  delays  and  reso¬ 
nance  positions  found  for  scattering  matrices  gen¬ 
erated  from  both  FE  and  variational  LHSF  are  listed 
in  table  1 ;  the  lifetime  matrix  eigenvalues  for  the  cur¬ 
rent  variational  LHSF  calculation  are  plotted  versus 
energy  in  fig.  1 .  The  main  difference  between  the  two 
calculations  js  the  appearance  of  two  weak  reso- 

"  In  a  previous  paper  [3]  we  used  Cl  to  denoie  ihe  vibrational 
angular  momentum  qua  mum  number,  in  analogs  lo  the  nota¬ 
tion  of  refs.  [16.17]  However,  we  have  also  used  Cl  [1.3]  to 
denoie  the  quantum  number  for  the  component  of  the  sys¬ 
tem's  toial  angular  momentum  along  the  direction  of  the  vec¬ 
tor  which  connects  the  center  of  mass  of  a  pair  of  ihe  system  's 
atoms  to  the  third  atom  (and  asymptotically  corresponds  to 
the  helicity  rotational  quantum  number).  Since  these  two  an¬ 
gular  momentum  components  are  in  general  distinct,  we  will 
for  clarity  use  the  symbol  K  in  this  paper  to  denote  the  first  one 
(i.e.  the  vibrational  angular  momentum  component),  while 
continuing  to  use  Cl  for  the  second. 
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Table  1 

Resonance  characteristics  for  PK.2  •’  potential  energy  surface 


J 

Assignment 

Current  results 

FEM 

lifetime  ( fs ) 

RPO" 

E  (eV) 

E  (eV) 

lifetime  t  fs  1 

E(eV) 

0.  1 

(0,  0°,0) 

0.61 

7 

061 

0.655 

1 

(0.  1‘.0) 

0.74 

6 

C,  1 

O 

►J 

O 

O 

0.85 

3 

0.847 

0.934 

0. 1 

(l.0°,0) 

0.97 

41 

0.971 

42 

0.975 

1 

(0,  3'.  0) 

0.97 

6 

0,  1 

(0,  4°,  0 ) 

1.07 

2 

1 

(t.  l'.O) 

1.08 

18 

0,  1 

( 1.  2°,  0) 

1.17 

7 

1.170 

1.175 

0.  1 

(2. 0",  0) 

1.38 

46 

1.382 

50 

t  .366 

1 

(2,  l'.O) 

1.47 

35 

0 

•> 

1.51 

7 

1.542 

0,  1 

(2.  2°,  0) 

1.56 

:o 

1.56 

1 

(2,3'. 0) 

1.65 

5 

•'  Ref.  [4 ], 

bl  Finite  element  results,  ref.  (3). 

‘ 1  Resonant  periodic  orbit  results,  ref.  [21], 


nances  at  energies  of  1.07  and  1.51  eV,  which  were 
not  previously  reported.  The  first  is  assigned  the  la¬ 
bel  (0.  4°,  0);  however,  the  resonance  at  1.51  eV,  in¬ 
dicated  by  an  unlabeled  arrow  in  fig.  1 ,  does  not  seem 
to  correspond  to  the  energy  of  an  expected  state  of 
metastable  linear  Hj  and  as  such  will  remain  unla¬ 
beled.  The  lifetimes  of  these  resonances  vary  greatly; 
the  long-lived  ones  at  0.969,  1.381,  and  1.56  eV  cor¬ 
respond  to  Feshbach  resonances  and  have  lifetimes 
of  41,  46,  and  20  fs,  respectively,  while  the  weaker 
peaks  correspond  to  shape  or  barrier  resonances  [21] 
and  have  an  average  lifetime  of  6  fs.  In  both  calcu¬ 
lations,  the  A,  and  E  symmetries  show  the  same  res¬ 
onance  energies  and  lifetimes,  with  more  numerical 
noise  present  in  the  E  calculation  due  to  the  larger 
number  of  states,  and  no  resonance  structure  is  found 
in  the  A 2  symmetry. 

The  J- 0  LSTH  surface  resonance  energy  and  life¬ 
times  from  the  current  calculations  and  the  previous 
FE  calculations  are  listed  in  table  2.  and  the  present 
lifetime  matrix  eigenvalues  are  displayed  in  fig.  2. 
The  assignment  of  these  states  is  the  same  as  for  those 
found  for  the  PK2  surface.  In  addition,  there  is  a  high 
energy  resonance  at  1.72  eV,  which  corresponds  to 
(2,  4°,  0);  however,  the  lifetime  analysis  near  this 
energy  is  obscured  '  •  numerical  noise,  so  this  energy 
is  less  reliable  than  the  other  resonance  energies. 
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The  energy  spacings  of  each  of  the  vz  =  0,  r:  =  2  and 
t’:  =  4  series  of,  resonances  suggest  that  a  resonance 
with  assignment  ( 1, 4°,  0)  should  exist,  for  the  LSTH 
surface,  at  the  position  indicated  in  fig.  2.  So  far,  this 
resonance  has  not  been  found.  By  analogy,  an  equiv¬ 
alent  resonance  should  exist  (but  is  not  found)  for 
PK_2,  as  indicated  in  fig.  1.  The  latter  cannot  cor¬ 
respond  to  the  unlabeled  resonance  at  1.51  eV,  be¬ 
cause  of  the  insufficiently  large  spacing  between  the 
(0,  4°,  0)  and  (2,  4°,  0)  resonances  for  the  LSTH 
surface. 

The  J=  1  partial  wave  includes  resonance  states 
with  K=  1  in  addition  to  AT=0.  Interestingly,  the  life¬ 
time  matrix  analysis  of  the  J- 1 ,  /7=  1 ,  A2  symmetry 
yields  the  same  resonance  energies  as  those  found  for 
the  J=  0,  A,  symmetry  for  both  the  PK2  and  LSTH 
surfaces,  with  one  exception;  we  therefore  interpret 
these  resonances  as  K=0  states.  The  exception  is  that 
there  is  no  visib'e  A;  resonance  at  1.51  eV  for  either 
surface,  but  it  is  possible  that  numerical  noise  inter¬ 
feres  with  its  detection.  The  absence  of  a  discemable 
energy  shift  due  to  the  increase  in  J  is  consistent  with 
the  approximate  rotational  constants  for  linear  Hj 
[23];  estimates  of  the  magnitude  of  the  shift  in  going 
from  7=0  to  J=  1  yields  a  value  of  about  0.002  eV, 
which  is  small  compared  with  the  accuracy  to  which 
we  have  determined  the  resonance  energies. 
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Fig.  1.  Resonant  lifetime  as  a  function  of  energy  for  the  A,  sym¬ 
metry  of  the /=0  partial  wave  of  H  +  H2  (PIC!  surface ).  The  ab- 
ctssa  £  is  the  total  energy  and  the  ordinate  represents  the  reso¬ 
nant  eigenvalue  of  the  collision  lifetime  matrix.  The  vertical 
arrows  on  the  upper  abcissa  denote  the  energies  of  the  H :(v.j) 
states.  The  length  of  these  arrows  decreases  as  v  spans  the  values 
0  through  3.  The  numbers  0,  5,  and  10  define  a  labeling  for  the 
value  of  j.  The  energy  grid  used  for  these  lifetime  calculations  was 
0  001  eV  m  the  neighborhood  of  the  two  strongest  resonances 
(l,0°,0)  and  (2, 0°,  0)  and  0.01  eV  elsewhere.  The  labeling  of 
the  resonances  at  the  top  of  the  panel  is  described  in  the  text.  The 
downward  pointing  unlabeled  vertical  arrow  at  1 .5 1  eV  indicates 
an  unassigned  resonance.  The  downward  pointing  arrow  labeled 
( 1 . 4°,  0 )  corresponds  to  a  resonance  expected  on  the  basis  of  en¬ 
ergy  spacings  ( see  text )  but  not  found  in  the  present  calculations. 

Additional  resonances  appear  in  the  J=  1 ,  77=  1 , 
A,  and  J=  1,  77=0,  A2  partial  waves.  For  the  PK2 
surface.,  they  occur  at  0.74,0.97,  1.08,  1.47,  and  1.65 
eV  (fig.  3).  The  assignments  and  lifetimes  are  given 
in  table  1.  The  r,  and  v2  assignments  are  done  on  the 
basis  of  the  energy  spacings,  and  the  K  assignment 
on  the  basis  of  the  restrictions  imposed  by  the  values 
of  J  and  v2  and  the  evenness  of  v2  [3,25].  The  cor¬ 
responding  values  for  the  LSTH  surface  are  listed  in 
table  2  and  displayed  in  fig  4;  they  have  energies  of 
0.77,  1.00,  1.09,  1.22,  1.45,  and  1.63  eV.  The  strong 
resonances  (1,  1 1 ,  0 ),  (2,  1 ',  0)  and  (2,  3 1 ,  0)  were 
found  previously  with  approximate  models  using  the 
LSTH  surface  [16,17,23,24];  these  results  are  also 
given  in  table  2  for  comparison.  The  remaining  peaks 


in  the  lifetime  analysis  are  weak  and  have  not  been 
reported  before.  No  resonances  were  found  for  the 
J=  1,  77=0,  A,  symmetry  on  either  surface.  All  res¬ 
onances  seen  in  the  J=  1  A,  and  A:  symmetries  with 
a  particular  parity  are  also  seen  in  the  J=  1  E  sym¬ 
metry  of  the  same  parity.  Again,  the  E  results  are  of 
lower  accuracy  because  of  the  larger  number  of  states. 

Calculations  for  collinear  triatomic  systems  have 
previously  given  strong  indications  [19]  that  the 
characteristics  of  resonance  spectra  are  closely  re¬ 
lated  to  the  geometry  of  potential  energy  surfaces  in 
the  strong  interaction  region  of  configuration  space, 
and  that  it  may  be  possible  to  infer  such  geometry 
from  experimentally  observed  resonance  spectra.  We 
will  now  try  to  obtain  such  relation  with  the  J=  0.  1 
resonance  at  hand.  Examination  of  the  energy  spac¬ 
ings  between  consecutive  resonances  in  each  of  the 
series  (0,  v2, 0),  ( 1,  v2,  0),  and  (2,  v2 ,  0)  shows 
them  to  be  nearly  constant  with  respect  to  v ,  and  v2 
and  having  global  averages  of  0.1 04  ±0.0 13  and 
0. 1 03  ±  0. 0 1 5  eV  for  the  LSTH  and  PK2  surfaces,  re¬ 
spectively.  This  correlates  very  well  with  the  spac¬ 
ings  of  0.1 1  and  0.12  eV  predicted  from  the  corre¬ 
sponding  bending  force  constants  [6,26]  and  a 
harmonic  model. 

Examining  the  series  ( v, ,  0°,  0 )  for  vx  =  0,  1.2  fur¬ 
nishes  consecutive  resonance  energy  differences  of 
0.33  and  0.38  eV  for  LSTH  and  0.36  and  0.41  eV  for 
PK2,  whereas  a  harmonic  model  based  on  the  sym¬ 
metric  stretch  force  constant  predicts  constant  spac¬ 
ing  of  0.26  eV  for  LSTH  and  0.47  eV  for  PK2.  Not 
surprisingly,  a  symmetric  stretch  static  model  does 
not  fit  the  resonance  spectra  well. 

The  energy  shift  between  the  (0.  0°,  0)  LSTH  and 
PK2  resonances  should  depend  in  part  on  the  dif¬ 
ference  of  0.029  eV  between  the  corresponding  sad¬ 
dle  point  energies  [4,6].  The  observed  downward 
shift  of  0.04  eV  can  be  totally  accounted  for  by  the 
difference  in  zero  point  bend  energies  and  saddle 
point  heights;  this  method  of  accounting  does  not 
seem  to  be  physically  reasonable,  since  the  difference 
in  the  symmetric  and/or  asymmetric  stretching 
characteristics  of  the  potential  energy  surfaces  should 
also  contribute  to  this  shift. 

We  conclude  that  bending  mode  force  constants  in 
the  saddle  point  region  of  this  system  can  easily  be 
obtained  from  the  corresponding  resonance  level 
spacings,  but  that  static  characteristics  of  the  sur- 
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Table  2 

Resonance  characteristics  for  LSTH  *'  potential  energy  surface 


J 

Assignment 

Current  results 

FE”' 

lifetime  (fs) 

RPO° 

E  (eV ) 

SCSA" 

£<eV) 

CEQB') 

E  (eV) 

CSM 

£<eV) 

E(eV) 

lifetime  (fs) 

£  (eV ) 

0.  1 

(C,  0°.0) 

0.65 

11 

0.65 

11 

1 

(0,  I'.O) 

0.77 

9 

0,  1 

(0,  2°,  0) 

0.88 

10 

0.880 

10 

0. 1 

(1.0°,  0) 
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0,  1 
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1.72" 
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1.734 

•' Refs.  (5.6).  01  Ref.  [22] 

cl  Resonant  periodic  orbit  results,  ref.  [23]. 

Small  curvature  semiclassical  adiabatic  results,  ref.  [24], 

*'  Colli  near  exact  quantum  with  adiabatic  bend  results,  ref.  [  1 7  J. 
fl  Coupled  state  results,  ref.  [16]. 

*’  This  resonance  energy  is  less  accurate  than  the  rest.  See  text  for  details. 


Fig.  2.  Same  as  Tig.  1  but  for  the  LSTH  surface;  a  constant  energy 
gnd  of  0.0 1  eV  was  used  throughout. 


Fig  3.  Resonant  lifetime  as  a  function  of  energy  for  the  A.  sym¬ 
metry  of  the/*  I.  /7=0  partial  wave  of  H  +  H2  (PK2  surface). 
The  abcissa  and  ordinate  are  as  given  in  fig  I;  the  energy  gnd  is 
0.0 1  eV  throughout.  The  downward  pointing  labeled  arrows  have 
similar  meanings  to  the  ones  in  figs.  I  and  2. 
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E/eV 

Fig.  A.  Same  as  fig.  3  but  for  the  LSTH  surface. 


Table  3 

Empirical  resonance  selection  Tiles 


J 

n 

K 

A, 

Aj 

E 

0 

0 

0 

yes 

no 

yes 

1 

0 

I 

no 

yes 

yes 

1 

i 

0 

no 

yes 

yes 

1 

i 

1 

vcs 

no 

yes 

faces  are  inadequate  to  understand  the  stretch  mo* 
non  features  of  the  resonance  spectra. 

The  presence  or  absence  of  resonances  in  each  of 
the  partial  waves  examined  is  summarized  in  table 
3.  An  empirical  selection  rule,  satisfied  by  the  results 
of  that  table,  is  that  resonances  are  present  in  the 
r=E  symmetry  for  all  allowed  values  of  K,  and  in 
r=A„  A:  when  the  quantity  (  —  1  )n*K  equals  Xr, 
where  Xr=  1  ( -  1 )  for  7~=  A,  (A:).  This  result  can 
be  derived  from  a  simple  model.  According  to  it,  no 
resonances  in  they.  17,  T  partial  wave  can  exist  if  the 
scattering  wavefunction  vanishes  identically  for  all 
configurations  of  the  system  in  the  vicinity  of  the 
saddle  point  for  which  the  distances  of  the  two  end 
atoms  to  the  central  atom  are  equal.  This  is  physi¬ 
cally  reasonable  since  we  expect  the  resonance  scat¬ 
tering  wavefunction  to  have  large  density  for  sym¬ 
metric  displacements  of  the  system  around  the  saddle 
point. 


Let  us  now  show  how  this  model  leads  to  the  em¬ 
pirical  rule  just  mentioned.  Let  ABC  be  a  linear  tn- 
atom  for  which  atoms  A  and  C  are  identical.  Con¬ 
sider  a  bent  configuration  of  this  system  in  which  the 
distances  AB  and  BC  are  equal.  Let  Xr  be  the  eigen¬ 
value  of  the  operator  which  permutes  A  and  C,  which 
is  equal  to  1  for  7~=  A,  and  to  —  1  for  /'=A:.  This 
permutation  is  equivalent  to  the  product  of  the  in¬ 
version  operator  and  a  rotation  by  n  around  the  sys¬ 
tem’s  principle  axis  of  inertia,  and  *herefore  appli¬ 
cation  of  these  operations  multiplies  the  triatom 
wavefunction  by  ( —  1  )n  and  ( —  1  )*,  respectively. 
Since  the  wavefunction  of  the  initial  symmetric  con¬ 
figuration  cannot  by  assumption  vanish  identically 
if  a  resonance  is  to  exist,  we  must  have 
Xr=  ( —  1  )n*K,  QED.  For  a  system  in  which  all  three 
atoms  are  identical,  neither  of  the  two  degenerate  E 
symmetry  wavefunctions  is  necessarily  even  or  odd 
with  respect  to  two-atom  permutations,  and  when 
one  of  these  wavefunctions  is  subjected  to  this  per¬ 
mutation  the  result  is  a  linear  combination  of  both; 
thus  the  E  symmetry  should  display  the  resonances 
found  in  both  the  A,  and  A:  symmetry  results  of  the 
same  parity,  as  we  have  indeed  observed. 

The  existence  of  resonances  seems  also  to  require 
the  presence  of  minima  in  adiabatic  curves  as  a  func¬ 
tion  of  an  appropriate  reaction  coordinate  [18- 
20,27-33].  In  the  particular  case  of  hyperspherical 
coordinates,  one  examines  LHSF  energies  including 
adiabatic  correction  terms  as  a  function  of  p;  these 
correction  terms  are  large  for  the  H  +  K2  system.  Plots 
of  LSHF  energies  for  the  J=  0, 1  partial  waves  of  each 
parity  and  symmetry  do  indeed  show  minima  (even 
without  corrections)  for  the  symmetry-parity  com¬ 
binations  which  support  resonances,  but  not  for  those 
combinations  which  have  shown  no  resonances.  A 
physical  interpretation  of  the  resonances  depends,  in 
addition  to  the  symmetry  arguments  given  above,  on 
an  explanation  as  to  why  these  particular  combina¬ 
tions  of  irreducible  representation  and  inversion 
parity  yield  adiabatic  energy  versus  p  curves  with 
minima. 


4.  Summary 

Application  of  the  variational  method  for  calcu¬ 
lation  of  LHSF  to  the  J=  l  partial  wave  of  H  +  H: 
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yields  results  of  sufficient  quality  for  lifetime  matrix 
analysis  to  give  accurate  resonance  energies  and  life¬ 
times.  The  dependence  of  resonance  energies  on  the 
bending  mode  characteristics  of  the  system’s  poten¬ 
tial  energy  surface  is  explained  by  a  simple  harmonic 
model,  and  the  existence  of  patterns  of  resonant  be¬ 
havior  in  terms  of  the  irreducible  representations  of 
the  P3  permutation  group  is  interpreted  using  a  sim¬ 
ple  physical  picture. 
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The  bound  rovibranonal  slates  of  the  upper  manifold  of  the  two  lowest  electronic  states  of  Hj  have  been  calculated  using 
variational  and  hyperspherical  coordinate  propagation  methods,  neglecting  in  both  the  coupling  between  those  electronic  states. 
Inclusion  of  the  effect  of  the  geometric  phase  induced  by  the  conical  intersection  between  those  manifolds  (sometimes  referred 
to  as  the  molecular  Aharonov-Bohm  effect )  is  shown  to  change  significantly  the  number,  the  energies  and  the  wavefunctions  of 
those  bound  rovibralional  states.  Quantum  numbers  are  defined  which  permit  a  physical  understanding  of  these  changes. 


1.  Introduction 

The  Rydberg  spectrum  of  the  H3  system  has  been  extensively  studied  by  Herzberg  and  coworkers  [I  ].  Of 
particular  interest  is  the  experimental  discovery  of  a  long-lived  metastable  state  [1-3].  On  the  theoretical  side, 
investigations  have  been  restricted  to  the  calculation  of  electronic  energies  for  a  few  nuclear  geometries  [4- 
7],  but  the  complete  electronic  potential  energy  surfaces,  necessary  to  investigate  the  rovibrational  structures 
of  the  spectrum  and  to  compute  accurately  the  lifetimes  of  the  excited  states,  are  available  only  for  the  ground 
and  the  first  electronically  excited  states  (DMBE  potential  [8] ).  In  the  equilateral  triangular  nuclear  config¬ 
uration,  these  two  electronic  states  are  degenerate  and  their  electronic  wavefunctions  belong  to  the  :E  rep¬ 
resentation  of  the  D3h  group.  Displacement  away  from  the  equilateral  triangular  geometry  lifts  this  degeneracy 
and  generates  a  conical  intersection  between  two  Jahn-Teller  sheets  [  9  ] .  Whereas  the  lower  sheet  is  responsible 
for  H  +  H2  reactive  scattering  below  about  3  eV  [10-17],  the  upper  one  supports  rovibrational  quasi-bound 
states,  which  can  predissociate  by  rovibronic  coupling  to  the  ground  electronic  state  [18,19]. 

In  this  Letter,  we  assume  that  the  upper  Jahn-Teller  sheet  is  decoupled  from  the  lower  one  and  therefore 
supports  bound  rovibrational  states.  We  compare  two  methods  of  computing  these  states  on  the  DMBE  excited 
potential  energy  surface.  One  is  the  variational  method  of  Tennyson  and  Sutcliffe  [20.2!  ]  (referred  to  as  TS 
method  in  this  paper)  The  other,  described  in  section  2,  is  a  hypersphencal propagation  method  which  uses 
modified  Whitten-Smuh  coordinates  ( 22.23  ]  and  derives  from  reactive  scattering  theory  [10-1 3.24  ] .  It  gen¬ 
eralizes  earlier  molecular  bound  state  calculations  limited  to  J- 0  [25.26].  We  show  in  section  3  that  the  hy- 
persphencal  method  is  very  appropriate  to  the  H,  system  because: 

-  It  allows  easy  inclusion  of  the  full  permutation  symmetries  of  the  three  identical  nuclei,  whereas  the  TS 
method  only  allows  inclusion  of  the  permutation  symmetries  of  two  identical  atoms. 

-  It  permits  inclusion  of  the  effect  of  the  conical  intersection  on  the  phase  of  the  nuclear  wavefunction  [27- 
29],  This  effect  results  from  the  sign  change  of  the  electronic  wavefunction  as  one  follows  a  closed  path  in  the 
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nuclear  configuration  space  around  the  line  along  which  the  two  :H'  electronic  states  conically  intersect.  It  cor¬ 
responds  to  a  particular  case  of  Berry's  geometric  phase  [30]  which  has  been  experimentally  observed  in  the 
Na3  system  [31  ].  Since  the  total  electronuclear  wavefunciion  is  continuous  and  single  valued,  there  has  to  be 
a  compensating  sign  change  in  the  nuclear  part  of  the  wavefunciion.  which  can  be  included  easily  in  the  hy- 
perspherical  method.  The  effect  of  the  conical  intersection  on  the  phase  of  the  nuclear  wavefunction  is  some¬ 
times  referred  10  as  the  molecular  Aharonov-Bohm  effect  [28.29.32],  but  we  will  use  the  simpler  name  "geo¬ 
metric  phase"  in  the  following. 


2.  Hyperspherical  method 


Let  A*,  A#,  Ar  be  the  atoms  of  the  system,  and  (A,  v ,  k)  be  any  cyclic  permutation  of  (a,  (3,  y ).  r,  is  the  mass- 
scaled  [33]  intemuclear  vector  for  the  diatom  A„A„  and  Rx  the  mass-scaled  vector  of  A^  with  respect  to  the 
center  of  mass  of  A^A^  The  hyperspherical  method  uses  the  hyperradius  p=  (Rl  +  ri) 1/2  to  describe  the  global 
size  of  the  triatomic  system  and  a  set  of  five  angles  £  to  describe  its  shape  and  orientation  in  space  [10- 
13,22,23,33.34].  In  this  paper,  we  will  neglect  all  spin-orbit  and  spin-spin  interactions.  In  the  Born-Oppen- 
heimer  approximation,  the  electronuclear  wavefunction  can  be  written  as  a  product  of  the  electronic  part 
which  we  choose  to  be  real,  and  the  nuclear  part.  The  latter  can  be  factored  into  a  nuclear  spin  part  and  a 
spacial  part  <pJsinr.  J  is  the  total  nuclear  angular  momentum  quantum  number.  Si  its  projection  onto  a  lab- 
oratory-fixed  axis,  /7the  parity  with  respect  to  the  inversion  of  nuclear  coordinates  and  Tthe  irreducible  rep¬ 
resentation  of  the  nuclear  permutation  group  (P3)  to  which  P Jx,nri  the  electronuclear  wavefunction  excluding 
the  nuclear  spin  part,  belongs: 

\pjMnr_ yj'.tnr^'  £)  ipt(qt\p,  £)  .  ( I  ) 

qc  refers  to  the  set  of  all.  spacial  and  spin,  electronic  coordinates.  ipJl{nr  is  an  eigenfunction  of  the  nuclear 
motion  Hamiltonian: 


1  >1  +  JL 

dpP  dp  lppl 


+  »■(/>.  o . 


(2) 


where  p  is  the  three-bodv  reduced  mass,  A  the  grand  canonical  angular  momentum  and  I'  the  Born-Oppen- 
heimcr  electronic  potential  energy  function.  The  nuclear  function  \pJ',nr  is  expanded  in  a  basis  of  local  hy¬ 
perspherical  surface  functions  (LHSF)  <pJn*nr\ 


vJ"nr(p,  £)  =  -Jr  I  £;  p )  . 

p  n 

The  LHSF  are  defined  as  the  eigenfunctions  of  the  fixed  hyperradius  nuclear  Hamiltonian: 
(— 7  +  Hp.  £)j  't>jnunr{*:p)=tjnnr(p)  <t>J”nr(Z.p)  . 


(3a) 


(3b) 


The  coefficients  F inr  in  eq.  ( 3a )  are  solutions  of  a  set  of  coupled  differential  equations  in  p.  w  hich  we  solve 
using  piece-wise  diabatic  bases  [10.34],  For  assumed  values  of  the  rovibrational  energies,  the  solutions  are 
propagated  forward  and  backward  from  small  and  large  p  values  where  they  have  negligible  amplitudes.  The 
energy  is  scanned  iteratively  until  the  quantization  condition  that  the  forward  and  backward  solutions  match 
smoothly  at  an  intermediate  value  of  p  is  reached. 

In  the  present  study,  we  use  the  Whitten-Smith  [22]  definition  of  the  five  angular  coordinates  £  as  modified 
by  Johnson  [23].  Three  Euler  angles  (aPy)  specify  the  orientation  of  the  body  frame  in  space.  The  axes  of 
this  frame  lie  along  the  principle  axes  of  inertia:  the  Z  axis  is  parallel  to  r,  x  R „  and  the  X  axis  is  associated  to 
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the  smallest  moment  of  inertia  and  is  oriented  such  that  rxx^  0.  Two  angles  (d,<px)  describe  the  shape  of  the 


molecular  triangle  and  are  defined  by 

r,.v  =^cos(n/4  — 0/2  )  sin(<p</2)  . 

(4a) 

rir=  — p  sin  (it/ 4  —  0/2 )  cos(^/2 )  , 

(4b) 

^=pcos(it/4-0/2)cos((Pi/2)  , 

(4c) 

Riy=psin(K/4-d/2)  sin(<px/2)  . 

(4d) 

The  ranges  for  these  angles  are  O^0^it/2  and  0^<px  <  27t.  0=0  corresponds  to  the  symmetric  top  configuration 
(an  equilateral  triangle  for  three  identical  particles)  in  which  the  principal  axes  of  inertia  X and  Y are  undefined. 
The  grand  canonical  angular  momentum  is  given  explicitly  by  [22,23] 


^J=-4AJ 


1  a 

sin  2030 


iin2eh 


4i ft  cos  0  *  _3_ 
+  sin20  2  <fyx 


2(P-PZ)  , 

cos20 


sin  0 
cos20 


(Ji+J  l). 


(5) 


where  3Z  is  the  body-fixed  Z  component  of  the  total  angular  momentum  J,  and  =Jx±'Jy. 

Eq.  (3b)  is  solved  variaiionallv  by  expansion  in  a  body-fixed  basis  xJJ££  built  with  products  of  simple  an¬ 
alytical  functions  [13]: 


X2X  =  exp(i n,9>i)/„(0)  Di,K(apy)  ,  (6) 

Djuk  is  a  Wigner  rotation  matrix  [35]  and  n¥  is  integer  or  half  of  an  odd  integer.  /„,(0)  are  simple  trigonometric 
functions,  such  ihat  the  LHSF  have  correct  behaviors  near  the  singularities  of  the  kinetic  energy  operator  0=0 
and  nl 2.  In  practice,  the /„  can  be  chosen  as  the  functions  cos (flg0)  or  sin (n<?0),  with  ng  integer  or  half  odd 
integer,  in  terms  of  which  the  hypersphencal  harmonics  ( whose  0  dependence  is  usually  written  as  a  polynomial 
in  cos0)  can  be  written  (eq.  (31 )  in  ref.  [36],  eqs.  ( 20 )— ( 23 )  in  ref.  [37]  or  eq.  (32)  in  ref.  [38] ). 

We  now  focus  attention  on  the  special  case  of  three  identical  nuclei  and  we  describe  how  to  build  electro- 
nuclear  wavefunctions  <pJMnr  which  are  bases  for  the  irreducible  representations  of  the  permutation  group  of 
the  nuclei  (P3).  The  operations  of  this  group  correspond  to  simple  changes  in  <px  (which  are  related  to  the 
isomorphism  between  P}  and  Clv)  as  indicated  in  table  1.  If  6*,  ( =  ±  1 )  is  the  symmetry  of  the  electronic 
wavefunction  with  respect  to  the  v+-k  permutation,  then  the  linear  combinations  defined  by 


,  .tn  ,r  /  ,  \J+K  +  2fvJM.-K  1 7 1 

r«l»,r  — —  1  )  Ant.- 1*1,1  • 

give  electronuclear  wavefunctions  {fJMnr  ( eq.  ( 1 ) )  with  the  ( =  ±  1 )  symmetry  with  respect  to  the  v*-k 
permutation. 

If  theie  is  no  conical  intersection  between  electronic  states,  the  electronic  wavefunction  <pt(qt;p,  O  belongs 
to  a  one-dimensional  representation  of  the  nuclear  permutation  group  ( At  for  =  +  1,  or  A2  for  t',K=  —  1 ). 
Table  2  indicates  how  the  total  angular  momentum,  the  parity  and  the  irreducible  representation  l~  of  P3  to 
which  'PJ'inr  belongs  determines  the  set  of  quantum  numbers  n^. 


Table  1 

Effect  of  permutations  of  the  nuclei  on  the  angle  ti 


Permuuuon 

mm 

p 

P  *  ’ 

r  r«i 

p  4) 

•  r« 

p  4) 
r  Iw 

P 

*  AM 

value  of  *,  ” 

»i  +  2n/3 

»4  +  4lt/3 

2* -9, 

2it/3-fi 

”  P.~  n  the  identity  permutation,  1,1  refers  to  the  cyclic  permutation  cl.  *’  P^.  refen  to  the  cyclic  permuuuon  /.vie— k/.v. 
* '  P„  refers  to  the  pair*  is«  permutation  of  nuclei  i  and  j. 

' 1  The  changes  in  are  true  modulo  2n.  since  Vj  must  remain  in  the  range  (0,  2it). 
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Table  2 

Choice  of  n,  for  each  parity  /7and  irreducible  representation  f of  the  nuclear  permutation  group  P, 


n 

r«i 

even,  v  uhout  phase  •’ 

A, /A, 

3m  al 

odd,  with  phase  *’ 

E 

3  mz  1  al 

even,  with  phase*1 

a,/a2 

3 m  +  ]  d> 

odd,  without  phase  *’ 

E 

3m±  j  <" 

*’  Without  consideration  of  the  geometric  phase  due  to  the  conical  intersection. 

*'  With  consideration  of  the  geometric  phase  due  to  the  conical  intersection. 

*’  r is  the  irreducible  representation  of  Pj  to  which  <pJMnr  (see  text  and  eq.  ( 1 ) )  belongs. 
4>  m  is  a  non-negative  integer. 


If  there  is  a  conical  intersection  between  electronic  states  for  equilateral  triangular  configurations  of  the  nu¬ 
clei  and  if  the  geometric  phase  is  taken  into  account,  one  can  show  [27-29]  that  in  the  vicinity  of  the  conical 
intersection  (6=0),  the  <px  dependence  of  the  Bom-Oppenheimer  electronic  wavefunction  is  given  by 


Wc^cos(<Px/2)  v/f 1  -sin(^/2)  wV 

(8a) 

or 

vrcaicos(9»J/2)  v/f,+sin(^/2)^f‘ 

(8b) 

where  (vrf1,  ^f!)  are  two  degenerate  /^-dependent  but  ^-independent  states  at  9=0  which  form  a  basis  for  the 
E  irreducible  representation  of  P3  being  symmetric  for  the  v~k  permutation  and  y/E2  antisymmetric). 
Although  permutations  of  the  nuclei  can  only  change  the  sign  of  y/t,  these  Bom-Oppenheimer  electronic  wave- 
functions  do  not  belong  to  a  one-dimensional  irreducible  representation  of  P3  and  are  discontinuous  in  the 
internal  configuration  space  [39]  in  the  plane  </>x  =  0.  However,  continuous  electronuclear  wavefunctions  which 
belong  to  irreducible  representations  of  P3  can  be  built  if  the  new  set  of  np  indicated  in  table  2  is  used  for  the 
nuclear  wavefunctions. 


3.  Results 

Fig.  1  illustrates  the  main  features  of  the  electronic  potential  in  the  internal  configuration  space  defined  in 


Fig.  I.  Plot  of  the  DMBE  excited  electronic  potential  r  in  the 
internal  configuration  space  defined  in  ref.  [39]  along  the  plane 
<*j='it/2, 3x/2  (i.e.  Z^O).  In  this  space,  the  coordinates  (p.  6.  <pk ) 
defined  in  the  text  correspond  to  spherical  polar  coordinates  with 
respect  to  the  Yt  axis  of  the  figure.  This  axis  is  also  the  one  along 
which  the  excited  DMBE  potential  conically  intersects  the  lower 
one.  The  equipotentials  are  equally  spaced  by  0.25  eV  in  the  range 
[3.  5  eV],  The  contours  for  V  =3  and  4  eV  are  specifically  indi¬ 
cated.  The  distances  on  the  A",  and  Yj  axes  are  in  bohr  Along 
constant  Yj  lines.  V  shows  the  usual  "V"-shaped  behaviour  char¬ 
acteristic  of  conical  intersections.  The  approximate  constancy  of 
the  spacing  between  the  equipotentials  in  this  figure  is  a  man¬ 
ifestation  of  this  linear  dependence.  Equipotentials  on  cuts  along 
other  planes  containing  the  Yj  axis  look,  in  the  vicinity  of  this 
axis,  very  similar  to  the  ones  displayed  in  this  figure,  i.e.  I  has  a 
local  nearly  cylindrical  symmetry  around  Yj. 
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ref.  [  T9  ] .  It  has  a  quasi-cvlindrical  symmetry  around  the  axis  of  that  figure  (6=0)  which  corresponds  to 
the  axis  of  the  conical  intersection  and  to  the  local  minima  on  the  fixed  p  spheres.  It  has  an  absolute  minimum 
for  p=  2.6  bohr.  0  =  0  corresponding  to  an  energy  of  2.72  eV  with  respect  to  the  bottom  of  the  ground  electronic 
state  H:  well.  In  the  vicinity  of  that  minimum,  the  potential  increases  steeply  and  almost  linearly  as  a  function 
of  6.  but  more  slowly  as  a  function  of  p. 

Table  3  compares  the  rovibrational  energy  levels  for  this  potential  obtained  by  the  hyperspherical  method 
and  the  TS  variational  method  without  consideration  of  the  geometric  phase. 

The  hyperspherical  method  uses  20  ne  values,  between  4  (A,  and  A2  symmetry)  and  8  (E  symmetry)  |n,| 
values  (eqs.  (6)  and  (7) ),  between  6  (A,  or  A2  symmetry)  and  12  (E  symmetry)  LHSF  (eq.  (3)).  The  LHSF 
have  been  computed  at  typically  50  p  values  between  1.5  and  6.5  bohr.  The  convergence  of  the  LHSF  and 
rovibrational  energies  is  of  the  order  of  10-4  eV.  The  compactness  of  the  hyperspherical  expansion  comes  from  • 
the  quasKylindrical  symmetry  of  the  potential  around  the  0=0  line  (small  number  of  n9  values)  and  from 
the  steep  increase  of  the  potential  as  a  function  of  0  (small  number  of  LHSF). 

The  TS  method  uses  a  body  frame  with  its  Z  axis  in  the  direction  of  R2  and  computes  the  bound  states  vari- 
ationally  by  expansion  on  a  product  basis  of  two  Morse-like  functions  (in  R2  and  q)  for  the  radial  part  and 
of  associated  Legendre  functions  for  the  angular  part.  The  optimized  parameters  of  the  Morse  potential  which 
we  chose  are  indicated  in  table  4.  Nearly  1400  such  product  functions  have  been  used  for  each  J ,  each  inversion 
parity  /7and  each  of  the  two  symmetries  for  the  permutation  of  the  two  identical  atoms  v  and  k.  This  unusually 
large  number  of  basis  functions  (only  880  such  functions  were  used  to  get  fully  converged  results  on  H3  in 
ref.  [40  J )  is  required  by  the  shape  of  the  potential  and  the  sudden  change  of  its  derivative  in  the  vicinity  of 
the  conical  intersection  axis.  Table  3  shows  that  the  convergence  of  the  energy  levels  is  always  worse  with  the 
TS  method  than  with  the  hyperspherical  method.  The  quality  of  the  TS  calculation  for  J=  1  odd  parity  is  not 


Table  3 

Bound  slate  energies  without  consideration  of  the  geometric  phase  *' 


mV  61 

J=0  c< 

J=  1  even  panty 

O 

J~  1  odd  parity 

C) 

000 

3.7210 

A, 

3.7218 

3.7283 

Aj 

3.7294 

3.7264 

E 

3.7276 

1  00 

3.9216 

A» 

3.9223 

3.9284 

Aj 

3.9297 

3.9266 

E 

3.9281 

200 

4.1067 

A| 

4.1073 

Aj 

4  1145 

4.1114 

E 

4.1131 

300 

4.2759 

A, 

4.2766 

4.2817 

A, 

4.2849 

4.2802 

E 

4.2831 

400 

4.4282 

A| 

4.4301 

4.4336 

Aj 

4.4386 

4.4322 

E 

4.4398 

500 

4.5621 

Aj 

4.5734 

4.5665 

A, 

4.5803 

4.5656 

E 

4.5894 

0  1  1 

4.2886 

E 

4.2886 

4.2955 

E 

4.2956 

4.2971 

A, 

4.2975 

4.2969 

A2 

4.2972 

4.2904 

£ 

4.2908 

1  1  1 

4.4533 

E 

4.4533 

4.4596 

E 

4.4598 

4.4610 

A, 

4.4618 

‘  4.4608 

Aj 

4.4615 

4.4550 

E 

4,4557 

2  1  1 

4.5980 

E 

4.5983 

4.6036 

E 

4  6048 

4.6049 

A| 

4.6083 

4.6047 

A2 

4.6093 

4.5996 

E 

4.6028 

3  1  1 

4.7212 

E 

4.7261 

E 

4.7349 

4.7272 

A| 

4.7370 

4.7270 

A2 

4.7355 

4.7225 

E 

020 

4  6806 

A, 

4  6813 

4.6871 

Aj 

4.6893 

4.6842 

E 

4.6878 

“  The  energy  is  in  eV  and  us  origin  corresponds  to  the  bottom  of  the  ground  electronic  state  of  the  isolated  H:  molecule. 

1,1  Quantum  numbers  used  to  classify  the  states  (see  text ) 

"  The  left  column  gives  the  hyperspherical  method  results  and  the  right  column  the  TS  method  results.  The  central  column  gives  the 
irreducible  representation  of  the  permutation  group  of  the  nuclei  to  which  the  special  pari  of  the  nuclear  wavefunction  belongs. 
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Table  4 

Optimized  parameters  of  the  Morse-like  functions  in  /?,  and 


O,  (au)  ■' 

ai,  ( au  1  " 

re  (au)  *’ 

J= 0 

0.230 

0.0130 

1.96  “* 

./  =  1 

0.262 

0.0100“' 

2.01  01 

J=  0 

0.262" 

0.01,22° 

2.09° 

Jss  1 

0.232° 

0.0102  ° 

2.32° 

These  parameters  are  defined  in  eqs.  (19)  and  (20)  of  ref.  (20). 
“  Parameters  for  the  Morse-like  functions  in  Rk. 

°  Parameters  for  the  Morse-like  functions  in  rA. 


as  good  as  the  TS  calculation  for  7=0  since  the  global  size  of  the  basis  has  been  kept  constant  instead  of  being 
doubled.  For  a  given  total  angular  momentum  and  parity,  the  quality  of  the  TS  results  decreases  as  the  energy 
increases,  and  in  particular,  states  diffuse  along/)  (corresponding  to  high  t'i  values,  see  below)  are  poorly  rep¬ 
resented.  This  suggests  that  different  sets  of  optimized  parameters  of  the  Morse-like  functions  should  be  used 
for  compact  and  diffuse  states. 

The  hyperspherical  method  can  be  compared  with  the  TS  method  from  computational  and  formal  points  of 
view: 

-  The  hypersphercal  method  requires  less  memory:  smaller  basis  sets  can  be  used  for  the  variational  solution 
of  the  two-dimensional  LHSF  equation  (see  eq.  (3b))  than  for  the  three-dimensional  variational  solution  of 
the  bound  states  in  the  TS  method.  However,  the  hyperspherical  method  required  about  two  times  more  CPU 
time  than  the  TS  method,  since  the  computation  of  the  LHSF  has  to  be  repeated  many  times,  but  did  not  exceed 
40  min  of  total  CPU  time  on  an  SCS-40  for  a  typical  run  7=0,  A,  plus  E  permutation  symmetries.  In  addition, 
the  hyperspherical  method  does  not  involve  adjustable  parameters  which  have  to  be  optimized  in  the  TS  method. 

-  The  bases  used  in  the  TS  method  to  expand  the  bound  state  wavefunciions  do  not  have  the  P3  permutation 
symmetry,  but  only  the  P2  symmetry  of  two  identical  nuclei.  As  a  result,  plots  of  the  bound  state  wavefunctions 
show  that,  even  in  the  7=0  case  where  the  energy  convergence  is  better  than  10~3  eV.  the  shape  of  the  TS 
wavefunctions  do  not  exhibit  the  correct  symmetry  properties  of  a  system  of  three  identical  particles,  whereas 
they  are  imbedded  in  the  LHSF  basis  used  in  the  hyperspherical  method.  Moreover,  the  TS  method  does  not 


> 


® 

LLI 


Fig.  2.  Rovibronic  energy  levels  associated  to  the  first  electroni¬ 
cally  excited  state  of  Hj.  The  full  lines  are  the  levels  including  the 
effect  of  the  geometric  phase  while  the  dashed  ones  exclude  that 
effect.  The  quantum  numbers  t, ,  t,  and  /  are  defined  in  the  text. 
The  origin  for  the  energy  scale  is  the  bottom  of  the  isolated  ground 
electronic  H-  potential  energy  curve.  These  levels  are  for  thei=0 
states,  but  the  y=  I  levels  are  nearly  degenerate  with  them,  the 
splitting  being  of  the  order  of  10~J  eV.  Their  nuclear  permuta¬ 
tion  symmetries  depend  on  J  and  on  ihe  parity  17.  as  well  as 
whether  the  geometric  phase  is  or  is  not  included  (see  tables  3 
and  5).  There  are  two  levels  for  each  of  the  sets  of  quantum  num¬ 
bers  (t>,  =  0.  rj  =  /=|)  and  ( r,  =  I .  t:  =  /  =  j ) .  which  would  be  de¬ 
generate  if  the  potential  were  exactly  cylindrical)  symmetric 
around  the  Y\  axis  (see  text  and  fig.  I  ). 
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Table  5 

Hypersp  lencal  method  energy  levels  including  the  effect  of  the  geometric  phase  *“’ 


l.r;/*’ 

J  =  0 

7=  1  even  parity 

J=  1  odd  parity 

o  !  ! 

4.0315  (F.) 

4.0286  (E) 

4.0256  (A,  I 

4.0243  (A,) 

4.0284  (E) 

1  !  ! 

4.2049  (E) 

.  4.2114(E) 

4.2087  (A,)  d’ 

4.2076  (Aj) 

4.2113  (E) 

■  2H 

4.3710  (E) 

4.3769  (E) 

4.3744  ( A, )  d> 

4.3734  (Aj) 

4.3768  (E) 

31! 

4.5189  (E) 

4.5243  (E) 

4.5220  (A,)  a> 

4.5210  (Aj) 

4.5241  (E) 

* ! ! 

4.6468  (E) 

4.6517  (E) 

4  6496  (A,) ai 

4.6487  (Aj) 

4.6515  (E) 

0  4  5 

4.5005  ( A , )  <* * 

4  5071  (Aj) 

4  5050  (E) 

4.5700  (Aj) 

4.5768  (  A, ) 

4.5733  (E) 

4.6425  (A,)  ■'  ■ 

4.6484  (Aj) 

4.6466  (E) 

4.7177  (Aj) 

4.7237  ( A, )  o’ 

4.7223  (E) 

* 1  The  energy  is  in  eV  and  ns  origin  corresponds  to  the  bottom  of  the  ground  electronic  state  of  the  isolated  H  j  molecule. 

The  irreducible  representations  are  the  ones  for  the  permutation  groupof  the  nuclei  to  which  yJ"nr belongs. 

' 1  Quantum  numbers  used  to  classify  the  stales  ( see  text ). 

dl  Levels  with  A,  symmetry  are  included  for  completeness,  but  are  forbidden  by  the  Pauli  principle. 


permit  inclusion  of  the  geometric  phase  due  to  the  conical  intersection  whereas  the  hypersphencal  method  does. 

Fig.  2  and  table  5  show  the  important  modifications  of  the  bound  rovibrational  energies  when  the  geometric 
phase  is  included  tn  the  hyperspherical  calculation.  These  changes  can  be  understood  if  one  defines  quantum 
numbers  to  the  bound  states  of  tables  3  and  5  by  modeling  the  nuclear  wavefunction  in  the  following  way 

-  We  retain  a  single  term  in  the  expansion  of  the  bound  states  in  the  LHSF  basis  (eq.  (3a) ).  This  Bom- 
Oppenheimer-type  approximation,  also  used  to  model  reactive  scattering  resonances  [24],  is  very  accurate  in 
the  present  case  where  the  frequency  associated  to  the  hyperspherical  mode  is  smaller  than  those  of  the  fixed- 
p  bending  modes:  the  resulting  bound  state  energies  are  shifted  by  less  than  0.4  meV.  This  approximation  sug¬ 
gests  that  we  define  the  quantum  number  v ,  associated  with  the  hyperradial  motion  as  the  number  of  nodes 
of  the  hyperradial  function  FJ,nr{p)  (eq.  (3a) ).  This  mode  corresponds  to  the  breathing  normal  mode  in  the 
limit  of  small  amplitude  vibrations,  but  in  the  present  case,  it  can  have  large  amplitudes  with  an  excitation  as 
large  as  r,  =  5  ( table  3  ). 

-  We  assume  that  the  fixed-p  bending  vibration  has  small  amplitude,  so  that  the  wavefunction  is  concentrated 
near  8=0.  This  approximation  is  reasonable  due  to  the  sleep  increase  of  the  potential  as  a  function  of  6.  It 
suggests  that  we  neglect  the  asymmetric  top  coupling  elements  in  the  kinetic  energy  ( last  term  of  eq.  ( 5  ) )  and 
the  dependence  in  the  potential.  The  (non-symmetrized)  LHSF  can  then  be  factored  as 

<P'"n  =  exp(in,$>J  gltl(d,p)  DJMK(a,P,y)  ,  (9a) 

where  gl},  is  defined  by 


•  I 


The  actual  energy  values  given  in  tables  3  and  5  are  calculated  accurately  ,  this  model  is  used  only  to  assign  quantum  numbers  to  these 
levels. 
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L  w>-\edd  d0 


j  &  ,/(£/>)  =  £«/*(*>)- 


/./(y-n)-i*2  ^  ik 

\  pp2 


— ;  g, 

w)  J 


(9b) 


Eq.  (9b)  is  the  small-0  limit  of  eq.  (3b)  (see  also  eq.  (5)).  /  quantizes  the  absolute  value  of  the  vibrational 
angular  momentum  in  a  new  body  frame,  which  is  an  Eckart  frame  associated  to  the  equilibrium  position  of 
the  nuclei  in  the  equilateral  triangular  configuration  [41  ],  and  is  given  by  l=\np—  J AT| .  v2  is  the  bending  vi¬ 
brational  quantum  number  and  is  defined  by  analogy  with  the  two-dimensional  harmonic  oscillator  such  that 
the  number  of  9  nodes  of  is  J(n2-/)  [42],  Vi  and  /  are  both  integers  when  the  geometric  phase  is  not- 
considered  and  become  both  half  odd  integers  when  it  is  taken  into  account.  If  the  pote  tial  were  a  harmonic 
function  of  6,  the  bound  state  energies  would  increase  linearly  with  tij  for  each  v{  value.  Although  the  potential 
is  an  approximate  linear  function  of  6,  tables  3  and  5  indicate  that  the  dependence  of  the  bound  state  energies 
on  t>2  is  not  far  from  linear.  Therefore,  as  shown  in  fig.  2,  each  of  the  levels  with  the  geometric  phase  (v2  half 
odd  integer)  is  almost  half  way  in  energy  between  two  consecutive  ones  without  this  phase  (i>2  integer). 

The  quantum  numbers  v2  and  l  defined  above  are  closely  related  to  the  ones  (n  and  j)  defined  in  ref.  [31  ] 
for  the  analysis  of  the  geometric  phase  in  the  22E'  Na3  excited  state.  However,  due  to  important  differences 
in  the  shapes  of  the  electronic  potentials  (minimum  for  equilateral  triangular  configurations  in  the  present 
excited  H3  state,  but  for  distorted  configurations  [9]  in  the  Na3  potential  used  in  ref.  [31  ] ),  the  dependence 
of  the  bound  state  energies  on  these  quantum  numbers  is  different  in  the  two  systems. 

Due  to  the  Pauli  principle  and  to  the  symmetries  of  the  nuclear  spin  wavefunction  with  respect  to  interchange 
of  the  identical  nuclei,  the  only  allowed  electronuclear  wavefunctions  'PJMnr  (eq.  ( 1 ) )  have  A2  or  E  nuclear 
permutation  symmetries,  and  they  correspond  to  quartet  and  doublet  nuclear  spins  respectively.  The  number 
of  such  levels  which  satisfy  the  Pauli  principle  and  their  spin  syn...ietries  change  significantly  when  the  effect 
of  the  geometric  phase  is  included. 


4.  Conclusions 

We  have  described  a  new  hyperspherical  propagation  method  for  the  calculation  of  bound  rovibrational  states. 
This  method  is  well  adapted  to  systems  of  three  identical  particles,  because  it  allows  easy  inclusion  of  the  full 
permutation  symmetries  of  the  system  and  of  the  effect  of  conical  intersections  on  the  phase  of  the  nuclear 
wavefunction. 

We  have  shown  that,  in  the  case  of  the  bound  rovibrational  states  in  the  first  electronically  excited  state  of 
Hj,  the  geometric  phase  results  in  bending  modes  having  half  odd  integer  quantum  numbers  and  in  important 
changes  of  the  rovibrational  state  energies  and  of  their  symmetry  properties.  In  the  following  paper  [43],  we 
study  the  influence  of  the  geometric  phase  on  the  chemical  reaction  which  occurs  in  the  ground  electronic  state 
of  H3. 
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The  effect  on  the  H  +  H  2  reaction  of  the  geometric  phase  induced  by  the  conical  intersection  between  the  two  lowest  electronic 
states  of  Hj  is  investigated  by  an  accurate  quantum  mechanical  computation  up  to  2.6  eV  above  the  bottom  of  the  ground  state 
H]  electronic  potential  well  for  the  total  angular  momentum  J^Q.  The  main  effects  of  the  inclusion  of  the  geometric  phase  are  a 
sign  change  in  the  reactive  scattering  matrix  and  modifications  in  the  nuclear  permutation  symmetries  of  the  calculated  reso¬ 
nances.  Cross-sections  to  which  non-reactive  processes  contribute  (ortho— ortho  and  para  — para)  are  significantly  modified, 
whereas  the  others  (ortho— para  and  para— ortho )  are  not. 


1.  Introduction 

Most  of  the  quantum  theoretical  studies  on  the 
H(2S)  +  H2(  )  reaction  have  used  the  Bom-Op- 

penheimer  approximation  (see  for  instance  refs.  ( I  — 
12] ),  and  assumed  that  the  reaction  occurs  on  the 
single  ground  electronic  potential  energy  surface.  This 
approximation  is  expected  to  be  quite  accurate  be¬ 
low  about  2.6  eV  of  total  energy  (with  respect  to  the 
bottom  of  the  H2(‘I*  )  potential  well)  since  this  is 
0. 1  eV  below  the  energy  of  the  minimum  of  the  first 
excited  potential. 

However,  a  complication  neglected  in  all  the  pre¬ 
vious  numerical  studies  [1-12]  arises  from  the  fact 
that  the  ground  electronic  state  conically  intersects 
the  first  excited  one  for  equilateral  triangular  con¬ 
figurations  of  the  nuclei.  There  is  a  sign  change  of  the 
electronic  wavefunction  as  one  follows  a  closed  path 
in  nuclear  configuration  space  around  the  line  of  the 
conical  intersection.  Since  the  total  electronuclear 
wavefunction  is  continuous  and  single-valued,  there 
has  to  be  a  compensating  sign  change  in  the  nuclear 
part  of  the  wavefunction  [13-16],  This  sign  change 
is  a  particular  case  of  Berry’s  geometric  phase  [17] 
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and  is  sometimes  referred  to  as  the  molecular  Aha- 
ronov-Bohm  effect  [14-16,18].  We  use  the  expres¬ 
sion  “geometric  phase”  in  the  rest  of  this  paper.  In 
the  preceding  paper  [19],  referred  to  hereafter  as  I, 
we  have  shown  that  this  geometric  phase  completely 
modifies  the  energy  spectrum  and  the  permutation 
symmetry  properties  of  the  quasi-bound  rovibra- 
tional  states  of  the  first  electronically  excited  state. 
In  this  Letter,  we  study  the  effect  of  this  geometric 
phase  on  the  chemical  reaction  which  occurs  on  the 
ground  electronic  potential  energy  surface. 

It  has  been  shown  formally  that,  if  the  condition 
that  the  wavefunction  is  zero  in  a  certain  region  of 
nuclear  configuration  space  separating  different  ar¬ 
rangement  channels  (see  section  3)  is  fulfilled,  the 
only  effect  of  the  geometric  phase  is  to  produce  an 
interference  between  reactive  and  non-reactive  scat¬ 
tering  amplitudes  which  is  exactly  the  opposite  to 
what  it  would  be  without  consideration  of  the  con¬ 
ical  intersection  [16].  This  condition  is  likely  to  be 
satisfied  for  the  low  collision  energies  considered  in 
the  earlier  quantum  studies  of  this  reaction  [1-3]. 
but  a  numerical  study  including  the  geometric  phase 
exactly  is  required  to  find  out  if  this  condition  re¬ 
mains  valid  at  the  higher  energies  reached  using  more 
recent  methods  [4-12], 

To  include  the  geometric  phase  in  the  calculation 
of  the  scattering  states  of  the  system,  two  distinct  ap- 

581 


0009-2614/90/$  03.50  *©  Elsevier  Science  Publishers  B.V.  (  North-Holland  ) 


95 


« 


Volume  166,  number  5,6 


CHEMICAL  PHYSICS  LETTERS 


9  March  1990 


proaches  can  be  followed,  both  leading  to  the  same 
final  result: 

-  One  can  use  real  electronic  and  nuclear  wave- 
functions  with  compensating  sign  changes  and  in¬ 
clude  the  geometric  phase  simply  by  enforcing  ap¬ 
propriate  boundary  conditions  and  permutation 
symmetries  on  the  nuclear  wavefunction. 

-  Alternatively,  one  can  add  extra  complex  phase 
factors  to  the  electronic  and  nuclear  parts  of  the 
wavefunction  to  enforce  the  continuity  and  single¬ 
valuedness  of  each  of  them.  These  extra  phases  add 
to  the  nuclear  Schrodinger  equation  a  term  formally 
similar  to  a  vector  potential  associated  to  a  delta- 
function  magnetic  field  located  on  the  conical  inter¬ 
section  line  [14,15],  Similarly  to  the  Aharonov- 
Bohm  effect  [18],  this  magnetic  field  modifies  the 
interference  pattern  between  semi-classical  trajec¬ 
tories  passing  on  opposite  sides  of  the  conical  inter¬ 
section  line  [14]. 

In  this  paper,  we  use  the  first  of  these  two  ap¬ 
proaches  to  compute  the  quantum-mechanical  three- 
dimensional  reactive  scattering  matrix  for  the  total 
angular  momentum  7=0  and  for  total  energies  be¬ 
low  2.6  eV.  We  restrict  the  use  of  the  second  ap¬ 
proach  to  the  semi-classical  interpretation  of  the  re¬ 
sults.  Section  2  describes  briefly  how  the  proper 
boundary  conditions  and  symmetry  properties  of  the 
nuclear  wavefunction  can  be  easily  incorporated  in 
a  hyperspherical  formalism  which  uses  modified  [20] 
Whitten-Smith  coordinates  [21]  to  describe  the 
strong  interaction  of  the  three  atoms  and  symme¬ 
trized  hyperspherical  coordinates  [4,22]  derived 
from  Delves'  coordinates  [23]  to  describe  the  initial 
and  final  H+H2  arrangements.  In  section  3,  we  de¬ 
scribe  in  detail  the  effect  of  the  geometric  phase  on 
the  scattering  matrix  and  cross-sections  for  distin¬ 
guishable  and  undistinguishable  nuclei. 


2.  Method 

We  outline  here  the  main  features  of  the  hyper¬ 
spherical  method,  already  described  in  I.  As  in  that 
paper,  we  neglect  spin-orbit  and  spin-spin  interac¬ 
tions  throughout.  In  the  Bohr-Oppenheimer  ap¬ 
proximation,  the  total  wavefunction  is  a  product  of 
electronic  and  nuclear  parts.  The  latter  is  a  product 
of  nuclear  spin  and  nuclear  spacial  wavefunctions. 
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The  electronuclear  wavefunction  <RJMnr  excluding  the 
nuclear  spin  part  (see  eq.  ( 1 )  in  I)  is  chosen  to  be¬ 
long  to  an  irreducible  representation  /"of  the  nuclear 
permutation  group  ( P3 )  of  H}.  It  is  also  labelled  by 
the  total  nuclear  angular  momentum  quantum  num¬ 
ber  J ,  its  component  M  along  a  laboratory-fixed  axis 
and  the  nuclear  parity  77.  In  the  absence  of  the  geo¬ 
metric  phase,  we  make  the  usual  assumption  that  the 
electronic  wavefunction  belongs  to  the  A|  irreduc¬ 
ible  representation  of  the  nuclear  permutation  group 
(P3).  In  this  case,  the  spacial  part  \pJMnr  of  the  nu¬ 
clear  wavefunction  belongs  to  the  same  irreducible 
representation  r  as  'PJMnr.  In  the  presence  of  the 
geometric  phase,  although  'PJ*nr  still  belongs  to  the 
irreducible  representation  Tof  P3,  the  spacial  part  of 
the  nuclear  wavefunction  and  the  electronic  wave- 
function  do  not. 

Let  Aa,  Ap,  A,  be  the  atoms  of  the  system,  and 
(A,  v,  k)  be  any  cyclic  permutation  of  (a,  p,  y).  rx  is 
the  mass-scaled  [23]  intemuclear  vector  for  the  di¬ 
atom  A  .A*  and  Rx  the  mass-scaled  vector  of  A„  with 
respect  to  the  center  of  mass  of  AA*-  In  the  hyper¬ 
spherical  method,  the  Hamiltonian  is  written  in  terms 
of  the  hyperradius  p-(Rl+rl)'n  which  parame¬ 
trizes  the  global  size  of  the  triatomic  system  and  of 
a  set  of  five  angles  C  which  describe  its  shape  and  ori¬ 
entation  in  space  [4,5,10-12,19-24]: 


id 


In  this  expression,  p  is  the  three-body  reduced  mass, 
A  the  grand-canonical  angular  momentum  and  V  the 
Bom-Oppenheimer  electronic  potential.  The  spacial 
nuclear  scattering  wavefunction  \pJKtnr  is  an  eigen¬ 
function  of  the  nuclear  Hamiltonian  (eq.  ( 1 ) )  and 
is  expanded  in  a  basis  of  local  hyperspherical  surface 
functions  (LHSF ),<PJ“nr: 


VJMnr(p,  0  =  -J7l  X  Fnnr(p )  <t>JnMnr{  C\p)  ,  ( 2 ) 

P  n 


which  are  eigenfunctions  of  the  fixed-hyperradius 
Hamiltonian: 

(2 ^  +  v^o)<t>iMnriC;p) 

=tinr(P)<Pi"nriQ\P).  (3) 

The  coefficients  FJmnr  in  eq.  (2)  are  solutions  of  a 
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set  of  ordinary  coupled  differential  equations  in  p. 
which  we  solve  using  piece-wise  diabatic  bases 
[4.25], 

A  solution  ofeq.  (3)  efficient  for  all  values  of  the 
hyperradius  requires  a  division  of  the  internal  con¬ 
figuration  space  into  two  regions: 

-  Region  I  corresponds  to  small  values  of  the  hy¬ 
perradius  for  which  the  three  atoms  interact  strongly. 
We  use  the  modified  [20]  Whiiten-Smith  [21  ]  def¬ 
inition  of  the  five  angular  coordinates  in  this  region: 
three  Euler  angles  define  the  orientation  of  the  frame 
of  the  principal  axes  of  inertia  in  space,  and  the  an¬ 
gles  (0,  (fi)  specify  the  shape  of  the  molecular  tri¬ 
angle.  The  LHSF  are  computed  numerically  by  ex¬ 
pansion  in  a  product  basis  of  simple  trigonometric 
functions  in  9  and  ^  {see  eqs.  (6)  and  (7)  in  I). 
Table  1  indicates  how  to  choose  the  functions  of  <pi 
to  obtain  electronuclear  wavefunctions  <PJMnr  with 
correct  permutation  symmetries,  with  and  without 
consideration  of  the  geometric  phase. 

-  Region  II  corresponds  to  large  values  of  the  hy¬ 
perradius  for  which  the  system  has  separated  into  an 
atom  and  a  diatom.  The  nuclear  wavefunction  is  now 
localized  in  the  electronic  potential  valleys  associ¬ 
ated  to  each  of  the  corresponding  asymptotic  ar¬ 
rangements  and  its  amplitude  on  the  plateaux  sep¬ 
arating  the  arrangement  channels  is  negligible.  This 
localization  makes  the  expansion  of  the  LHSF  on  the 
delocalized  product  basis  used  in  region  I  inefficient 
and  suggests  the  use  of  symmetrized  hyperspherical 
coordinates  [4,22]  instead.  Therefore,  {0,  <px)  are 
replaced  in  region  II  by  {colt  y,,)  defined  by  (ox= 
2arctan (rA//?A)  and  yi=arccos(Rl-ri/Riri),  each 
in  the  range  0  to  it.  The  LHSF  are  now  expanded  in 
a  product  basis  [4]  of  Legendre  polynomials  in  cosy* 


and  vibrational-type  functions  in  a>x.  Since  product 
bases  associated  to  different  arrangements  do  not 
overlap  in  region  II,  the  LHSF  which  include  the 
geometric  phase  differ  from  the  ones  which  exclude 
it  only  by  simple  changes  in  the  signs  of  the  pieces 
of  the  wavefunction  within  each  arrangement  chan¬ 
nel.  The  geometric  phase  can  be  included  straight¬ 
forwardly  in  region  II  since  it  does  not  change  the 
overlap  and  interaction  matrices. 

The  LSTH  Bom-Oppenheimer  electronic  poten¬ 
tial  energy  surface  has  been  used  [26].  This  surface 
was  chosen,  rather  than  the  DMBE  [27]  one  used 
for  the  bound  state  calculations  of  paper  I,  because 
we  already  had  accurate  scattering  results  [4]  for  it 
which  served  to  test  the  validity  of  the  new  hyper¬ 
spherical  method  described  here.  Recent  results  [  7,9  ] 
indicate  that  the  use  of  the  DMBE  potential  does  not 
significantly  change  the  final  scattering  matrix  re¬ 
sults.  The  boundary  between  region  I  and  region  II 
was  set  at  p— 6  bohr.  Surface  functions  were  com¬ 
puted  at  20  values  of  p  between  2  and  6  bohr  in  re¬ 
gion  I  and  3 1  values  of  p  between  6  and  1 2  bohr  in 
region  II.  The  results  shown  below  have  been  ob¬ 
tained  with  1 1 56  product  functions  ( 34  for  each  of 
the  two  angular  coordinates  9  and  <px )  in  region  I  and 
156  product  functions  in  region  II.  We  verified  that 
the  convergence  of  the  scattering  matrix  elements 
with  respect  to  the  size  of  the  product  basis  was  of 
the  order  of  1%  by  comparing  with  a  smaller  cal¬ 
culation  involving  900  product  functions  in  region 
I.  65  LHSF  for  A,  and  A2  permutation  symmetries 
and  1 30  LHSF  for  E  symmetries,  with  energies  ac¬ 
curate  to  within  approximation  10-J  eV,  have  been 
used  in  the  expansion  of  the  wavefunction.  The  uni- 
tarity  of  the  resulting  scattering  matrix  was  always 


Table  I 

Basis  in  for  expansion  of  the  7=0  LHSF 


A," 

A," 

E” 

without  phase  bl 

cos  (3n0,J 

sin  ( 3/i0i ) 

cos[(3/i±  1  )0i] 

with  phase c) 

cos[(3/t  +  ])0il 

sin[  (3n  + j)0i] 

cos[(3/t+  j)0i) 

*’  Irreducible  representation  of  the  permutation  group  of  the  nuclei  to  which  the  electronuclear  wavefunction  9*,*",rb  jongs.  We  choose 
the  component  of  the  E  irreducible  representation  which  is  symmetric  with  respect  to  the  A,— A,  permutation  (tee  text),  n  is  a  non¬ 
negative  integer. 

*’  Excluding  consideration  of  the  geometric  phase.  In  this  case,  the  electronic  wavefunction  is  assumed  to  belong  to  the  A,  irreducible 
representation  of  P,. 

*’  Including  consideration  of  the  geometric  phase. 
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better  than  1%.  except  for  the  E  symmetry  above  2.4 
eV  for  a  few  columns  and  rows  associated  to  highly 
excited  reactants  or  products. 

3.  Results 

We  first  discuss  the  effect  of  the  geometric  phase 
on  the  scattering  matrices  associated  to  electrons 
dear  wavefunctions  'pJ"nr  which  are  bases  for  ir¬ 
reducible  representations  of  the  permutation  group 
of  the  nuclei.  Then,  we  switch  to  the  scattering  ma¬ 
trices  assodated  to  distinguishable  particles.  This 
representation  affords  a  simple  understanding  of  the 
effect  of  the  geometric  phase  on  the  H  +  H:  reaction. 

Figs.  1  and  2  show  some  transition  probabilities 
for  the  total  angular  momentum  7=0  and  for  elec- 
tronuclear  wavefunctions  belonging  to  irre¬ 

ducible  representations  of  Pj.  Comparison  of  both 
rows  of  figs.  1  and  2  shows  that  the  inclusion  of  the 
geometric  phase  induces  important  changes  in  the 
qualitative  features  of  the  transition  probabilities  as¬ 
sociated  to  irreducible  representations  of  the  per¬ 
mutation  group,  except  for  the  transitions  from  para 
to  ortho  (and  ortho  to  para  because  of  the  micro¬ 
scopic  reversibility)  hydrogen  in  the  E  irreducible 
representation  which  remain  almost  unchanged  (last 
column  of  fig.  2). 

One  interesting  modification  due  to  the  inclusion 
of  the  geometric  phase  is  the  exchange  of  qualitative 
features  between  the  transition  probabilities  associ¬ 
ated  to  the  A i  and  A2  irreducible  representations  (fig. 
1 ).  In  particular,  the  resonances,  which  correspond 
to  A  i  wavefunctions  without  geometric  phase,  ap¬ 
pear  on  the  Aj  probabilities  when  this  phase  is  in¬ 
cluded.  This  exchange  results  from  a  similar  one  in 
the  nodal  structure  of  the  nuclear  wavefunctions.  In¬ 
deed,  the  last  of  ref.  [  4  ]  shows  that  resonance  wave- 
functions  are  non-zero  for  near  linear  configurations 
of  Hj  where  two  intemuclear  distances  are  equal. 
These  configurations  correspond  to  the  half  planes 
^=it/3,  it,  5x/3,  which,  according  to  table  1,  are 
nodal  planes  for  the  A2  wavefunctions  in  the  absence 
of  the  geometric  phase  and  for  the  A,  wavefunctions 
when  this  phase  is  included.  Therefore,  the  reso¬ 
nances  must  correspond  to  A,  electronuclear  wave- 
functions  '¥JMnr  without  the  geometric  phase  and  to 
A2  wavefunctions  with  this  phase  for  7=0. 
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Fig.  1.  Transition  probabilities  as  a  function  of  energy  for  y=  0 
without  (upper  row)  and  with  (lower  row)  consideration  of  the 
geometric  phase.  The  left  column  corresponds  to  an  electronu¬ 
clear  wavefunction  <pJt,nr  which  belongs  to  the  A,  irreducible 
representation  of  the  permutation  group  of  the  nuclei  and  to 
the  transition  H  +  Hj(u»0,>»0,  m/»0)-H-t-Hj<v'  =  l..;'=0, 
mrm0).  The  right  column  corresponds  to  an  electronuclear 
wavefunction  <pJMnr  which  belongs  to  the  A2  irreducible  repre¬ 
sentation  of  the  permutation  group  of  the  nuclei  and  to  the  tran¬ 
sition  H+Hj(uw0,>«  I,  m,=0)—H+Hj(u'  *1,/  *1,  m,  ®0). 
The  lower  abscissa  is  the  total  energy  and  the  upper  abscissa  the 
reagent  relative  translational  energy.  The  vertical  arrows  on  the 
upper  abscissa  denote  the  energies  of  the  Hj(i>,;=0)  states  and 
are  labelled  by  the  values  of  v.  The  squares  on  the  probability 
curves  indicate  the  points  for  which  the  scattering  calculations 
were  made. 

More  generally,  the  inclusion  of  the  geometric 
phase  modifies  the  result  of  the  symmetry  analysis 
of  resonances  for  low  total  angular  momentum  de¬ 
scribed  in  the  last  of  ref.  [4].  This  paper  shows  that 
the  exchange  of  the  two  extreme  hydrogen  atoms  of 
the  near  linear  configurations  considered  in  the  pre¬ 
vious  paragraph  multiplies  the  spacial  part  of  the  nu¬ 
clear  wavefunction  by  a  factor  ( —  1  )n**,  where  K  is 
the  vibrational  angular  momentum  of  the  resonant 
state  and  17  its  parity.  If  Xr  is  the  eigenvalue  of  the 
operator  which  permutes  the  two  extreme  hydrogen 
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Fig.  2.  Transition  probabilities  as  a  function  of  energy  for  J= 0 
without  (upper  row)  and  with  (lower  row)  consideration  of 
the  geometric  phase,  for  the  transitions  H  +  Hj(t/=0,;  =  0, 
m,  =  0)  —  H  +  Hj(i-'  =«  1,/  =»0.  m,  =0)  (left  column)  and  H  + 
HilrwO.v^O,  m;=0)  — H  +  Hj(u'  =  I,/  =  1,  m,  =0)  (right  col¬ 
umn)  with  the  electronudear  wavefunction  <pJMnr  belonging  to 
the  E  irreducible  representation  of  P,.  For  the  other  details,  see 
the  caption  of  fig.  1. 

atoms  in  the  electronudear  wavefunction  'PJMnr  ( + 1 
for  the  A,  irreducible  representation,  - 1  for  A2  and 
±  1  for  the  E  doubly  degenerate  irreducible  repre¬ 
sentation,  according  to  which  component  is  consid¬ 
ered),  the  condition  for  having  a  resonance  without 
consideration  of  the  geometric  phase  is 
Xr=  ( -  1  )n*K,  since,  in  this  case,  we  assume  that 
the  electronic  part  of  the  wavefunction  is  symmetric 
for  the  permutation  of  the  hydrogen  nuclei.  If  now 
the  geometric  phase  is  included,  the  condition  for 
having  a  resonance  becomes  Xr=  (  - 1  )r7+K* 1  since 
the  electronic  wavefunction  is  now  antisymmetric 
with  respect  to  the  nuclear  permutation.  Therefore, 
the  symmetry  assignments  of  the  resonances  are  ex¬ 
changed  for  the  A,  and  A2  irreducible  representa¬ 
tions,  but  remain  unchanged  for  the  E  irreducible 
representation,  which  has  doubly  degenerate  sym¬ 
metric  and  antisymmetric  components. 


We  now  assume  that  the  three  hydrogen  atoms  are 
distinguishable  particles.  Fig.  3  shows  that  the  non¬ 
reactive  (direct)  and  reactive  (exchange)  transition 
probabilities  are  almost  not  affected  by  the  inclusion 
of  the  geometric  phase.  Only  slight  changes  (less  than 
10%)  appear  in  the  probabilities  above  2eV  for  low 
excited  states  of  the  reactants  and  products  (these 
changes  however  can  become  larger  than  10%  for 
small  transitions  between  rotationally  excited  reac¬ 
tant  and  product  states).  The  strong  effect  observed 
in  fig.  1  in  the  transition  probabilities  associated  to 
the  irreducible  representations  of  the  permutation 
group  of  the  nuclei  results  mainly  in  a  change  of  n  in 
the  phase  of  the  reactive  scattering  matrix  elements , 
whereas  the  phases  of  the  non-reactive  matrix  ele¬ 
ments  and  the  norms  of  the  reactive  and  non-reac¬ 
tive  scattering  matrix  elements  are  only  slightly 
modified  by  the  geometric  phase.  These  numerical 
results  validate  the  conclusion  of  ref.  [16]  and  the 
assumption  on  which  it  rests  over  a  wide  energy 
range.  Indeed,  ref.  [16]  shows  formally  that,  if  the 
wavefunction  is  zero  in  the  vicinity  of  the  half  plane 
tpi  =  it,  then  the  only  effect  of  the  geometric  phase  is 
to  change  the  sign  of  the  reactive  scattering  matrix. 
This  result  should  be  also  valid  for  non-zero  total  an¬ 
gular  momenta. 

Ref.  [14]  suggests  a  semiclassical  picture  for  the 
effect  of  the  geometric  phase  on  the  reactive  prob¬ 
abilities,  in  terms  of  the  modification  of  the  inter¬ 
ference  pattern  between  trajectories  passing  on  op¬ 
posite  sides  of  the  conical  intersection  axis  (0=0). 
However,  the  quantum  results  suggest  that,  since  only 
the  phases  of  the  reactive  scattering  matrix  elements 
are  significantly  changed,  almost  all  trajectories 
should  pass  on  the  same  side  of  the  conical  inter¬ 
section,  namely  the  one  nearest  the  minimum  energy 
path.  We  performed  quasi-classical  trajectory  cal¬ 
culations  for  J=0  with  a  sample  of  400  trajectories 
for  each  collision  energy.  Below  2.6  eV.  we  found  in¬ 
deed  that  for  the  ground  rovibrational  initial  state 
(v-j=0)  only  one  trajectory  (at  2.4  eV)  passes  on 
the  side  of  the  conical  intersection  opposite  to  the 
minimum  energy  path.  However,  rotational  excita¬ 
tion  of  the  reactants  can  increase  this  number.  10% 
of  the  reactive  trajectories  pass  on  the  side  opposite 
to  the  minimum  energy  path  at  2.5  eV  for ;=  10  and 
15.  This  suggests  that,  as  observed  in  the  quantum 
results,  the  effect  of  the  geometric  phase  can  become 
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Fig.  3.  Non-reactive  (left)  and  reactive  (right)  transition  probabilities  for  H  +  Hj(r=0,;=0,  m,  =  0)  —  H  +  H2(t/  =  I,/  =0.  m,  =0). 
The  dots  and  the  continuous  lines  refer  to  the  calculations  excluding  and  including  the  geometric  phase  respectively.  For  the  other  details, 
see  the  caption  of  fig.  1 . 


more  important  for  some  transitions  between  rota- 
tionally  excited  reactants  and  products.  Finally,  the 
number  of  trajectories  passing  on  the  side  of  the  con¬ 
ical  intersection  opposite  to  the  minimum  energy 
path  increases  strongly  above  2.8  eV  and  reaches  40% 
of  the  total  number  of  reactive  trajectories  for  j— 0 
at  4  eV.  This  indicates  that  the  geometric  phase 
modifies  both  the  norms  and  the  phases  of  the  scat¬ 
tering  matrix  elements  in  this  high  energy  range. 
However,  inclusion  of  the  electronic  coupling  to  the 
first  electronically  excited  state  becomes  necessary  to 
obtain  quantitatively  correct  results  in  this  energy 
range. 

Consideration  of  the  Pauli-antisymmetrized  (de¬ 
fined  in  the  second  ref.  [  1  ] )  7=0  partial  cross-sec¬ 
tions  allows  us  to  estimate  qualitatively  how  the 
Pauli-antisymmetrized  integral  cross-sections  are 
modified  by  inclusion  of  the  geometric  phase.  The 
previous  discussion  indicates  that  the  para  to  ortho 
or  ortho  to  para  antisymmetrized  cross-sections  are 
almost  not  modified  for  energies  below  2.6  eV  by  in¬ 
clusion  of  the  geometric  phase,  since  only  the  reac¬ 
tive  scattering  amplitude  appears  in  the  expressions 
of  the  antisymmetrized  cross-sections  for  these  tran¬ 
sitions  (seeeqs.  (5.39)  and  (5.40)  in  the  second  ref. 
[  1  ] ).  Therefore,  this  phase  cannot  be  the  reason  for 
the  important  discrepancies  between  experimental 
and  theoretical  integral  cross-sections  below  1.4  eV 
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for  para  to  ortho  transitions:  the  strong  resonant  pat¬ 
terns  which  are  present  on  the  experimental  integral 
cross-sections  [28]  will  not  appear  on  the  corre¬ 
sponding  theoretical  results  [6,7,12]  even  if  the  geo¬ 
metric  phase  is  included.  However,  fig.  4  suggests  that 
antisymmetrized  cross-sections  can  be  significantly 
modified  by  the  introduction  of  the  geometric  phase 
for  para  to  para  or  ortho  to  ortho  transitions,  since 
the  way  reactive  and  non-reactive  scattering  ampli¬ 
tudes  interfere  is,  with  a  good  accuracy  below  2.6  eV, 
changed  to  its  opposite  (see  eqs.  (5.39)  and  (5.40) 
in  the  second  ref.  [  I  ] ).  This  change  is  not  very  im¬ 
portant  at  low  energy  when  the  reaction  probabilities 
are  small  compared  to  the  non-reactive  ones  (see  the 
000—  020  transition  probability  below  0.6  eV  in  fig. 
4),  but  it  can  modify  even  the  qualitative  features 
and  the  order  of  magnitude  of  the  antisymmetrized 
cross-sections  when,  at  higher  energy,  reactive  and 
non-reactive  scattering  amplitudes  become  of  the 
same  order  of  magnitude. 

The  effect  of  the  geometric  phase  on  the  details  of 
the  angular  distributions  of  the  ortho— ortho  and 
para— para  cross-sections  should  be  even  more  pro¬ 
nounced  than  on  the  corresponding  integral  cross- 
sections.  The  differential  cross-sections  for  these 
transitions  show  an  oscillatory  dependence  on  scat¬ 
tering  angle  [  1  ]  which  should  become  more  intense 
as  energy  increases.  These  oscillations  are  due  to  in- 
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Fig.  4.  Antisymmetrized  J= 0  partial  integral  cross-sections  for  the  transition  H  +  H2  (i  =  0.;  =  0,  m,  =  0)  —  H  +  H2(r'  =0. /  =  2,  m,  =0) 
(left)  and  H  +  H2  (i  =  0.7  =  0.  m,  =  0)  —  H  +  H2(i/’  =  1,/  =0.  m,  =0)  (right).  The  dots  and  the  continuous  lines  refer  to  the  calculations 
excluding  and  including  the  geometric  phase  respectively.  For  the  other  details,  see  the  caption  of  fig.  I . 


terferences  between  direct  and  exchange  scatterings, 
and  a  change  in  the  sign  of  the  exchange  scattering 
amplitude  should  make  constructive  interferences 
destructive  and  vice  versa  (16). 

4.  Conclusions 

The  present  numerical  study  of  the  geometric  phase 
in  the  H  +  H2  validates  the  conclusion  of  ref.  [  1 6  ]  in 
a  wide  energy  range  below  2.6  eV:  quite  accurate 
cross-sections  can  be  obtained  by  neglecting  the  geo¬ 
metric  phase  in  the  computations  of  the  reactive  and 
non-reactive  scattering  matrix  elements  and  by  in¬ 
cluding  it  a  posteriori  by  changing  the  signs  of  the 
reactive  scattering  matrix  elements.  This  sign  change 
can  modify  significantly  the  spin-averaged  cross-sec¬ 
tions  when  the  energy  is  high  enough  for  reactive 
transition  probabilities  to  be  non-negligible  com¬ 
pared  to  the  non-reactive  ones.  Quasi-classical  tra¬ 
jectory  calculations  indicate  that  this  approximate 
treatment  of  the  geometric  phase  becomes  inaccur¬ 
ate  at  the  higher  energies  for  which  the  Born-Op- 
penheimer  approximation  is  expected  to  break  down. 
Finally,  our  quantum  study  also  indicates  that  in¬ 
clusion  of  the  geometric  phase  changes  the  symmetry 
assignments  of  the  resonances. 
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We  have  performed  accurate  three-dimensional  quantum  mechanical  reactive  scattering  calculations  for  the  H  +  H:  system  on 
the  Caltech/JPL  Mark  fllfp  64  processor  hvpercube,  using  the  method  of  symmetrized  hyperspherical  coordinates  and  local 
hyperspherical  surface  functions.  The  results  and  timing  obtained  demonstrate  that  such  distributed  memory  parallel  architec¬ 
tures  are  competitive  with  the  CRAY  X-MP,  CRAY  2  and  CRAY  Y-MP  supercomputers  and  should  allow  the  study  of  larger, 
more  complicatedchemical  systems.  In  addition,  we  show  that  a  selection  rule  for  scattering  resonances  developed  previously  and 
tested  for7=0.  I  resonances  is  also  satisfied  by  the  7=2  resonances  obtained  in  the  present  calculations. 


I.  Introduction 

There  is  considerable  current  interest  in  perform¬ 
ing  accurate  quantum  mechanical  three-dimensional 
reactive  scattering  cross  section  calculations.  Accu¬ 
rate  solutions  have  until  recently  proved  to  be  dif¬ 
ficult  and  computationally  expensive  to  obtain,  in 
large  part  due  to  the  lack  of  sufficiently  powerful 
computers  [1-7].  Prior  to  the  advent  of  supercom¬ 
puters,  one  could  only  solve  the  equations  of  motion 
for  model  systems  or  for  sufficiently  light  atom-di¬ 
atom  systems  at  low  energy  [  1-4  ].  As  a  result  of  the 
current  development  of  efficient  methodologies  and 
increased  access  to  supercomputers,  there  has  been 
a  remarkable  surge  of  activity  in  this  field  [8-19], 
The  use  of  symmetrized  hyperspherical  coordinates 
[20]  and  of  the  local  hyperspherical  surface  function 
formalism  [8,9,21  ]  has  proven  to  be  a  successful  ap¬ 
proach  to  solve  the  three-dimensional  Schrodinger 
equation  [8,9,15,16],  However,  even  for  modest  re¬ 
active  scattering  calculations  the  memory  and  CPU 
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the  Ph.D.  degree  in  Chemistry  at  the  California  Institute  of 
Technology. 
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demands  are  so  great  that  CRAY-type  supercom¬ 
puters  will  soon  be  limiting  progress. 

Although  there  has  been  a  steady  improvement  in 
the  necessary  technologies  of  the  basic  logic  speeds 
of  computers,  there  is  little  prospect  of  substantially 
faster  single  processor  designs  in  the  near  future. 
Concurrent  supercomputers  are  a  natural  next  step 
in  meeting  the  need  for  both  increased  memory  and 
faster  CPU.  Individual  processors,  although  slower 
than  a  single  sequential  supercomputer  processor,  can 
be  connected  together  in  sufficient  number  to  make 
a  powerful  supercomputer.  Such  architectures  offer 
the  potential  to  obtain  large  increases  in  computing 
speed  by  simply  increasing  the  number  of  proces¬ 
sors.  The  actual  speed-up  depends  on  the  nature  of 
the  algorithm,  the  characteristics  of  the  processors, 
and  the  particular  way  these  communicate  with  each 
other.  The  algorithms  used  and  the  codes  developed 
on  sequential  machines  should  be  replaced  by  codes 
optimized  for  parallel  machines. 

The  essential  property  a  calculation  must  have  to 
be  efficiently  done  on  a  highly  parallel  computer  is 
that  it  be  decomposable  in  such  a  way  that  in  per¬ 
forming  it  almost  all  processors  should  be  computing 
efficiently  almost  all  of  the  time,  and  that  the  com¬ 
munication  time  between  the  processors  should  rep¬ 
resent  a  small  fraction  of  the  computation  time.  In 
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the  present  paper  we  show  how  quantum  mechanical 
reactive  scattering  calculations  can  be  structured  so 
as  to  fulfil  these  criteria. 

The  hvpercube  architecture  is  a  leading  design  for 
MIMD-type  (multiple  instruction  multiple  data) 
distributed  memory  parallel  architectures  based  on 
message  passing.  The  first  such  machine  was  devel¬ 
oped  by  Seitz  [22]  and  used  by  Fox  [23,24],  both 
at  Caltech.  We  have  created  efficient  codes  to  solve 
the  quantum  mechanical  equation  of  motion  for  re¬ 
active  collisions  of  an  atom  w;th  a  diatomic  mole¬ 
cule  using  a  hypercube  computer  of  this  type.  Very 
similar  codes  should  be  appropriate  for  other  MIMD 
distributed  memory  parallel  architectures. 

In  this  paper,  we  present  a  concurrent  algorithm 
for  calculating  local  hvperspherical  surface  functions 
(LHSF)  and  use  a  parallelized  version  [25]  of  John¬ 
son’s  logarithmic  derivative  method  [26],  modified 
to  include  the  improvements  suggested  by  Manolo- 
poulos  [27],  for  integrating  the  resulting  coupled 
channel  reactive  scattering  equations.  We  review  the 
formalism  briefly  in  section  2.  In  section  3  we  dis¬ 
cuss  the  parallel  algorithms  and  in  section  4  we  com¬ 
pare  the  results  of  scattering  calculations  on  the  Cal- 
tech/JPL  Mark  Illfp  64  processor  hypercube  for  the 
H  +  H2  system  7=0,  1,2  partial  waves  on  the  LSTH 
[28,29]  potential  energy  surface  with  those  of  cal¬ 
culations  done  on  a  CRAY  X-MP/48  and  a  CRAY- 
2.  Both  accuracy  and  performance  are  discussed,  and 
speed  estimates  are  made  for  the  Mark  Illfp  1 28  pro¬ 
cessor  hypercube  soon  to  become  available  and  the 
San  Diego  Supercomputer  Center  CRAY  Y-MP/864 
machine  which  has  just  been  put  into  operation.  We 
summarize  the  conclusions  in  section  5. 


2.  Methodology 

The  detailed  formulation  of  reactive  scattering 
based  on  hvperspherical  coordinates  and  local  vari¬ 
ational  hvperspherical  surface  functions  (LHSF)  is 
discussed  elsewhere  [8.9,1 5  ].  We  present  a  very  brief 
review  to  facilitate  the  explanation  of  the  parallel 
algorithms. 

For  a  triatomic  system,  we  label  the  three  atoms 
A„,  A9  and  Ar  Let  (/.,  v,  k)  be  any  cyclic  permuta¬ 
tion  of  the  indices  (a,  p,  y).  We  define  the  A  coor¬ 
dinates,  the  mass-scaled  [30]  intemuclear  vector  rx 
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from  A„  to  A«,  and  the  mass-scaled  position  vector 
Rx  of  \x  with  respect  to  the  center  of  mass  of  the  A,  A, 
diatom.  The  symmetrized  hyperspherical  coordi¬ 
nates  [20]  are  the  hyperradius  p=  (/?;  +r; ) l/:,  and 
a  set  of  five  angles  u>x,  yx,  6X,  0X  and  wx,  denoted  col¬ 
lectively  as  £i.  The  first  two  of  these  are  in  the  range 
0  to  rc  and  are  respectively  2  arctan  (rx/Rx)  and  the 
angle  between  Rx  and  rx.  The  angles  8X,  <t>x  arc  the  po¬ 
lar  angles  of  Rx  in  a  space-fixed  frame  and  to*  is  the 
tumbling  angle  of  the  Rx,  rx  half-plane  around  iis  edge 
Rx.  The  hamiltonian  Hx  is  the  sum  of  a  radial  kinetic 
energy  operator  term  in  p,  and  the  surface  Hamil¬ 
tonian  dx,  which  contains  all  differential  operators  in 
Ci  and  the  electronically  adiabatic  potential  V(p,  tox, 
yx).  (tx  depends  on  p  parametrically  and  is  therefore 
the  "frozen”  hyperradius  part  of  Hx. 

The  scattering  wave  function  '¥JMnr  is  labelled  by 
the  total  angular  momentum  7,  its  projection  vr  on 
the  laboratory-fixed  Z  axis,  the  inversion  paruy  77 
with  respect  to  the  center  of  mass  of  the  system  and 
the  irreducible  representation  r  of  the  permutation 
group  of  the  system  (P3  for  H  +  H2)  to  which  the 
electronuclear  wave  function,  excluding  the  nuclear 
spin  part  [31,32],  belongs.  It  can  be  expanded  in 
terms  of  the  LHSF  <pJMnry  defined  below,  and  cal¬ 
culated  at  the  values  pq  of  p: 

W"nr(p,  Ci>»  X  bi?r(p-pq)<t>i""r{^pq)  .  ( 1 ) 

It 

The  index  i  is  introduced  to  permit  consideration  of 
a  set  of  many  linearly  independent  solutions  of  the 
Schrodinger  equation  corresponding  to  distinct  ini¬ 
tial  conditions  which  are  needed  to  obtain  the  ap¬ 
propriate  scattering  matrices. 

The  LHSF  0iw/7r(Ci',  Pq)  and  associated  energies 
tinr(pq)  are  respectively  the  eigenfunctions  and  ei¬ 
genvalues  of  the  surface  Hamiltoniat-'ii.  The;,  are 
obtained  using  a  variational  approacn  [15]  The 
variational  basis  set  consists  of  products  of  Wigner 
rotation  matrices  DJ„a(0x,  0X.  tpx),  associated  Le¬ 
gendre  functions  of  yx  and  functions  of  cox  which  de¬ 
pend  parametrically  on  p„  and  are  obtained  from  the 
numerical  solution  of  one-dimensional  eigenvalue- 
eigenfunction  differential  equations  in  a>x  involving 
a  potential  related  to  V((j,  a >x,  yx). 

The  variational  method  leads  to  an  eigenvalue 
problem  with  coefficient  and  overlap  matrices 
hjnr(pq)  and  sjnr(pq)  and  whose  elements  are  five- 
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dimensional  integrals  involving  the  variational  basis 
functions. 

The  coefficients  bJ„?r(p\  pq)  defined  by  eq.  (1  ) 
satisfy  a  coupled  set  of  second-order  differential 
equations  involving  an  interaction  matrix  Jjnl  (p\ 
p„)  whose  elements  are  defined  by 

[Jjnr(p\Pq)  jr  =  I  V(p,  wx,  yj 

-(PJp)1  y(P^u)l,yi)\<t>iMnr(^pq)'>  .  (2) 

The  configuration  space  p ,  £,  is  divided  in  a  set  of  Q 
hyperspherical  shells  i  (4=1.  2,  ....  Q) 

within  each  of  which  we  choose  a  value  p„  used  in 
expansion  ( 1 ). 

When  changing  from  the  LHSF  set  at/?,  to  the  one 
at  /?,„ ,  neither  ,pJl'-,nr  nor  its  derivative  with  respect 
to  p  should  change.  This  imposes  continuity  condi¬ 
tions  on  the  bJ’!r  and  their  p  derivatives  at  /?=/?,+  ,, 
involving  the  overlap  matrix  Cjnr(pq+I,  p,)  be¬ 
tween  the  LHSF  evaluated  at  p,  and  p,+  , 

[  Cjnr(pq+ 1 .  p,)  ]!' 

=  <*J:,nrc*.Pn)\<pti,nr(^pq)>  .  (3) 

The  five-dimensional  integrals  required  to  evalu¬ 
ate  the  elements  of  hjnr,  sjnr,  Jjnr  and  Cjnr  are 
performed  analytically  over  6X  and  tpx  and  by  two- 
dimensional  numerical  quadratures  over  yx  and  cu^. 
These  quadratures  account  for  90%  of  the  total  time 
needed  to  calculate  the  LHSF  <pJ“nr  and  the  ma¬ 
trices  Jjnr  and  Cjnr. 

The  system  of  second-order  ordinary  differential 
equations  in  the  bJ£r  is  integrated  as  an  initial  value 
problem  from  small  values  of  p  to  large  values  using 
Manolopoulos’  logarithmic  derivative  propagator 
[27],  Matrix  inversions  account  for  more  than  90% 
of  the  time  used  by  this  propagator.  All  aspects  of  the 
physics  can  be  extracted  from  the  solutions  at  large 
p  by  a  constant  p  projection  [8,9,33]. 


3.  Parallel  algorithm 

The  computer  used  for  this  work  is  a  64  processor 
Mark  lllfp  hypercube.  Each  node  consists  of  two  in¬ 
dependent  Motorola  68020  microprocessors,  one  for 
computation  and  one  for  I/O,  and  four  megabytes  of 
dynamic  memory.  The  computation  microprocessor 
has  a  Motorola  68882  floating-point  arithmetic  co¬ 


processor  and  128  kilobytes  of  static  private  mem- 
ory.The  I/O  microprocessor  has  64  kilobytes  of  static 
private  memory.  An  additional  daughter  board  with 
a  pipe-lined  32-bit  floating  point  unit  based  on  the 
Weitek  XL  senes  of  chips  is  attached  io  each  node 
and  has  a  nominal  peak  speed  of  16  Mflops.  The 
crystalline  operating  system  (CrOS)-channel-ad- 
dressed  synchronous  communication  provides  the  li¬ 
brary  routines  to  handle  communications  between 
nodes  [24,34,35].  Program  development  is  done  on 
a  Motorok  68020-based  Counterpoint  workstation 
that  runs  on  UNIX.  The  programs  are  written  in  C 
programming  language  except  for  the  time-consum¬ 
ing  two-dii.iensional  quadratures  and  matrix  inver¬ 
sions,  .ich  are  optimized  in  assembly  language. 

The  hypercube  is  configured  as  a  two-dimensional 
array  of  processors.  The  mapping  is  done  using  bi¬ 
nary  Gray  codes  [24.36]  which  gives  the  Cartesian 
coordinates  in  processor  space  and  communication 
channel  tags  for  a  processor's  nearest  neighbors.  With 
a  distributed-memory  machine  like  the  hypercube, 
the  elements  of  a  large  matrix  of  data  must  be  dis¬ 
tributed  across  the  memory  of  all  the  processors.  This 
makes  it  possible  to  fully  utilize  the  large  memory 
available  and  facilitates  the  load-balancing  task  of 
keeping  most  of  the  processors  busy  doing  useful 
arithmetic  most  of  the  time.  The  parallelization  of 
scientific  codes  is  frequently  based  on  a  large  grain 
size  decomposition  of  the  task.  A  method  of  distrib¬ 
uting  the  global  matrix  among  the  processors  is  the 
first  choice  that  must  be  made  and  it  is  closeiy  re¬ 
lated  to  the  parallel  algorithm  chosen. 

We  mapped  the  matrices  into  processor  space  by 
local  decomposition.  Let  Nr  and  Nc  be  the  number  of 
processors  in  the  rows  and  columns  of  the  hypercube 
configuration,  respectively.  Element  A  ( /,  j)  of  an 
M x  A 1  matrix  is  placed  in  processor  row  P,  =  int  ( iSJ 
M)  and  column  where  int.v  means 

the  integer  part  of  x. 

The  parallel  code  implemented  on  the  hypercube 
consists  of  five  major  steps.  Step  one  constructs,  for 
each  value  of  /?,,  a  primitive  basis  set  composed  of 
the  product  of  Wigner  rotation  matrices,  associated 
Legendre  functions,  and  the  numerical  one-dimen¬ 
sional  functions  in  u>x  mentioned  in  section  2  and 
obtained  by  solving  the  corresponding  one-dimen¬ 
sional  eigenvalue-eigenvector  differential  equation 
using  a  finite  difference  method.  This  requires  tha- 
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a  subset  of  the  eigenvalues  and  eigenvectors  of  a  tri- 
diagonal  matrix  be  found. 

A  bisection  method  [37]  which  accomplishes  the 
eigenvalue  computation  using  the  TR1DIB  routine 
from  EISP.ACK.  [38]  was  ported  to  the  Mark  Illfp. 
This  implementation  of  the  bisection  method  allows 
computation  of  any  number  of  consecutive  eigen¬ 
values  specified  by  their  indices.  Eigenvectors  are 
obtained  using  the  EISPACK.  inverse  iteration  rou-' 
tine  TINVIT  with  modified  Gram-Schmidt  ortho- 
gonalization.  Each  processor  solves  independent  tri¬ 
diagonal  eigcnproblems  since  the  number  of 
eigenvalues  desired  from  each  tridiagonal  system  is 
small  but  there  are  a  large  number  of  distinct  tridi¬ 
agonal  systems.  To  achieve  load  balancing,  we  dis¬ 
tributed  subsets  of  the  primitive  functions  among  the 
processors  in  such  a  way  that  no  processor  computes 
greater  than  one  eigenvalue  and  eigenvector  more 
than  any  other.  These  large  grain  tasks  are  most  eas¬ 
ily  implemented  on  MIMD  machines:  S1MD  (single 
instruction  multiple  data)  machines  would  require 
more  extensive  modifications  and  would  be  less  ef¬ 
ficient  because  of  the  sequential  nature  of  effective 
eigenvalue  iteration  procedures.  The  one-dimen¬ 
sional  bases  obtained  are  then  broadcast  to  all  the 
other  nodes. 

In  step  two  a  large  number  of  two-dimensional 
quadratures  involving  the  primitive  basis  functions 
which  are  needed  for  the  variational  procedure  are 
evaluated.  These  quadratures  are  highly  parallel  pro¬ 
cedures  requiring  no  communication  overhead  once 
each  processor  has  the  necessary  subset  of  functions. 
Each  processor  calculates  a  subset  of  integrals 
independently. 

Step  three  assembles  these  integrals  into  the  real 
symmetric  dense  matrices  sjnr(pq)  and  hjnr{pq) 
which  are  distributed  over  processor  space.  The  en¬ 
tire  spectrum  of  eigenvalues  and  eigenvectors  for  the 
associated  variational  problem  is  sought.  With  the 
parallel  implementation  of  the  Householder  method 
[  39  ] .  this  generalized  eigensystem  is  tridiagonalized 
and  the  resulting  single  tridiagonal  matrix  is  solved 
in  each  processor  completely  with  the  QR  algorithm 
[40].  The  QR  implementation  is  purely  sequential 
since  each  processor  obtains  the  entire  solution  to 
the  eigensystem.  However,  only  different  subsets  of 
the  solution  '■re  kept  in  different  processors  for  the 
evaluation  of  .he  interaction  and  overlap  matrices  in 
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step  four.  This  part  of  the  algorithm  is  not  time-con¬ 
suming  and  the  straightforward  sequential  approach 
was  chosen.  It  has  the  further  effect  that  the  resulting 
solutions  are  fully  distributed,  so  no  communication 
is  required. 

Step  four  evaluates  the  two-dimensional  quadra¬ 
tures  needed  for  the  interaction  Jjnr(p\  pq)  and 
overlap  Cjnr(pq+,  ;  p„)  matrices.  The  same  type  of 
algorithms  are  used  as  were  used  in  step  two.  By  far, 
the  most  expensive  part  of  the  sequential  version  of 
the  surface  function  calculation  is  the  calculation  of 
the  large  number  of  two-dimensional  numerical  in¬ 
tegrals  required  by  steps  two  and  four.  These  steps 
are,  however,  highly  parallel  and  well  suited  for  the 
hypercube. 

Step  five  uses  Manolopoulos’  [27]  algorithm  to 
integrate  the  coupled  linear  ordinary  differential 
equations.  The  parallel  implementation  of  this  al- 
goriaim  is  discussed  elsewhere  [25].  The  algorithm 
is  dominated  by  parallel  Gauss-Jordan  matrix  in¬ 
version  and  is  I/O  intensive,  requiring  the  input  of 
one  interaction  matrix  per  integration  step.  To  re¬ 
duce  the  I/O  overhead  a  second  source  of  parallel¬ 
ism  is  exploited.  The  entire  interaction  matrix  (at  all 
p)  and  overlap  matrix  (at  all pq)  data  sets  are  loaded 
across  the  processors  and  many  collision  energies  are 
calculated  simultaneously.  This  strategy  works  be¬ 
cause  the  same  set  of  data  is  used  for  each  collision 
energy  and  because  enough  main  memory  is  avail¬ 
able.  Calculation  of  scattering  matrices  from  the  fi¬ 
nal  logarithmic  derivative  matrices  is  not  compu¬ 
tationally  intensive,  and  is  done  sequentially. 

The  program  steps  were  all  run  on  the  Weitek  co¬ 
processor  which  only  supports  32-bit  arithmetic.  Ex¬ 
perimentation  has  shown  that  this  precision  is  suf¬ 
ficient  for  the  work  reported  below.  The  64-bit 
arithmetic  hardware  needed  for  larger  calculations 
was  installed  after  the  present  calculations  were 
completed. 

4.  Results  and  discussion 

Accuracy.  Calculations  were  performed  for  the 
H  +  H2  system  on  the  LSTH  surface  [28,29]  for  par¬ 
tial  waves  with  total  angular  momentum  J- 0.  1.  2 
and  energies  up  to  1 .6  eV.  Flux  is  conserved  to  better 
than  1%  for  J=  0,  2.3%  for  J=  1  and  3.6%  for  J=  2 
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Fig.  1.  Probabilities  (a)  and  probability  differences  (b)  as  a 
function  of  total  energy  £  (lower  abscissa)  and  initial  relative 
translational  energy  Em  (upper  abscissa)  for  the  J=0  (0,  0. 
0)-.  (0, 0, 0)  A,  symmetry  transition  in  H  +  H2  collisions  on  the 
LSTH  potential  energy  surface.  The  symbol  ( v,  j.  Q )  labels  an 
asymptotic  state  of  the  H  +  H  2  system  in  which  v.j,  and  fl  are  the 
quantum  numbers  of  the  initial  or  final  H2  stales.  The  vertical 
arrows  on  the  upper  abscissa  denote  the  energies  at  which  the 
corresponding  H2(v,  j )  states  open  up.  The  length  of  those  ar¬ 
rows  decreases  as  v  spans  the  values  0,  I  and  2,  and  the  numbers 
0,  3,  and  10  associated  with  the  arrows  define  a  labelling  for  the 
value  of  j.  (a)  The  results  from  the  Mark  Hlfp  hypercube;  (b) 
differences  between  these  and  those  from  the  CRAY  X-MP/48. 
The  number  of  LHSF  used  was  36  and  the  number  of  primitives 
used  to  calculate  these  surface  functions  was  80. 


Fig.  2.  Same  as  for  fig.  1  except  for/=  1.  A,,  odd  parity  (77=  ! ). 

( 0,  0,  0 )  —  ( 0,  0.  2 )  transitions.  The  number  of  LHSF  used  was 
74  and  the  number  of  primitives  used  to  calculate  these  surface 
functions  was  1 52. 

for  ali  open  channels  over  the  entire  energy  range 
considered. 

To  illustrate  the  accuracy  of  the  32-bit  arithmetic 
calculations,  the  scattering  results  from  the  Mark  IUfp 
with  64  processors  are  shown  in  figs.  1.  2.  and  3  for 
7=0,  1,  2,  respectively,  in  which  some  transition 
probabilities  as  a  function  of  the  total  collision  en¬ 
ergy  E  are  plotted.  Also  shown  are  the  differences 
between  these  results  and  those  obtained  using  a 
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E/eV 

Fig.  3.  Same  as  for  fig.  1  except  for 7=2,  A,,  odd  parity  (77=  1 ), 
(0,  2,  I )  —  (0, 2,  I )  transition.  The  differences  plotted  in  (b)  are 
between  the  Mark  lllfp  hypercube  and  the  CRAY-2  results.  The 
number  of  LHSF  used  was  63  and  the  number  of  primitives  used 
to  calculate  these  surface  functions  was  1 36. 

CRAY  X-MP/48  and  a  CRAY-2.  These  differences 
do  not  excede  0.004  in  absolute  value  over  the  en¬ 
ergy  range  investigated.  The  effect  of  the  geometric 
phase  associated  with  the  conical  intersection  be¬ 
tween  the  two  lowest  electronic  potential  energy  sur¬ 
face  of  H3  [32]  is  not  included  in  these  results.  Much 
of  the  structure  in  the  transition  probability  curves 
is  due  to  the  underlying  resonances  [1,9,  16]  and 
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are  discussed  below.  The  two  sets  of  data  in  each  fig¬ 
ure  are  virtually  indistinguishable  on  the  scale  of  the 
plots. 

Analysis  of  7=  2  resonances.  Table  1  contains  a  list 
of  the  J=  2  resonance  energies  detected  from  the 
maxima  in  the  lifetime  versus  energy  curves,  calcu¬ 
lated  as  described  previously  [9,16],  as  well  as  their 
quantum  number  assignments,  permutation  and  in¬ 
version  symmetry  labels,  and  lifetimes.  The  per¬ 
mutation  symmetries  are  given  with  and  without  the 
inclusion  of  the  effect  of  the  geometrical  phase  (GP ) 
associated  with  the  conical  intersection  between  the 
two  lowest  electronic  state  potential  energy  surfaces 
[31,32],  The  energy  of  these  resonances  is  consis¬ 
tent  with  the  physical  model  for  the  selection  rule 
previously  developed  [16]  and  tested  with  the  7=0. 
1  resonances.  The  results  of  table  1  adds  additional 
credence  to  the  generality  of  that  rule.  According  to 
it,  if  GP  effects  are  ignored,  a  necessary  (but  not  suf¬ 
ficient  )  condition  for  resonances  to  occur  in  A ,  ( A 2 ) 
partial  waves  is  that  ( —  1  )n+K  be  equal  to  1  (  —  1 ). 
where  K  is  the  vibrational  angular  momentum  quan¬ 
tum  number,  whereas  they  are  permitted  in  E  partial 
waves  for  all  K.  To  include  the  GP  effect,  it  suffices 
to  interchange  A,  and  A2  in  this  selection  rule  [32], 
In  agreement  with  this  picture,  not  all  higher  energy 
7=2  resonances  which  are  allowed  by  this  rule  were 
detected. 

Timing  and  parallel  efficiency.  In  tables  2  and  3  we 
present  the  timing  data  on  the  64  processor  Mark 
IHfp,  a  CRAY  X-MP/48  and  a  CRAY-2,  for  both 
the  surface  function  code  ( including  calculation  of 
the  overlap  Cjnr  and  interaction  Jjnr  matrices )  and 
the  logarithmic  derivative  propagation  code.  For  the 
surface  function  code,  the  speeds  on  the  first  two 
machines  is  about  the  same.  The  CRAY-2  is  1.43 
times  faster  than  the  Mark  lllfp  and  1.51  times  faster 
than  the  CRAY  X-MP/48  for  this  code.  The  reason 
is  that  this  program  is  dominated  by  matrix-vector 
multiplications  which  are  done  in  optimized  assem¬ 
bly  code  in  all  three  machines.  For  this  particular  op¬ 
eration  the  CRAY-2  is  2.03  times  faster  than  the 
CRAY  X-MP/48  whereas  for  more  memory -inten¬ 
sive  operations  the  CRAY-2  is  slower  than  the  CRAY 
X-MP/48  [41  ].  A  slightly  larger  primitive  basis  set 
is  required  on  the  Mark  lllfp  in  order  to  obtain  sur¬ 
face  function  energies  of  an  accuracy  equivalent  to 
that  obtained  with  the  CRAY  machines.  This  is  due 
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Table  I 

/=2  resonances  and  their  characteristics  for  H3 


Assignment 

Permutation  sy  mmetry  *' 

Inversion 

symmetry 

fl” 

£(eV) 

Lifetime  ( fs ) 

without  GP  61 

with  GP 01 

(0. 0°.0) 

A„E 

A,.  E 

0 

0.65 

11 

(0.  I'.O) 

A-.  A„E 

AJt  A2,  E 

0.  1 

0.77 

9 

(0, 2°,  0) 

A|.  E 

A 2,  E 

0 

0.88 

10 

(0. 2,.0) 

A,.A„E 

A2,  A|,  E 

0,  1 

0.90 

10 

( 1, 0®.  0) 

A,.E 

a2,  e 

0 

0.98 

28 

(0,3',  0) 

A2,  A(,  E 

A|.  A2,  E 

0.  1 

1. 00 

8 

(1.  I'.O) 

A2»  A|,  E 

A|»  A2,  E 

0.  1 

1.09 

29 

(0.4°,0) 

A,.E 

A  2,  E 

0 

1.10 

5 

(0, 4L0) 

A„Aj,E 

a2,  A  j,  e 

0,  1 

1.12 

5 

( 1. 2J,  0) 

A,.A„E 

Aj.  A,.E 

0.  1 

1.22 

8 

(O.S'.O) 

A2,  A,,  E 

A|,  Aj,  E 

0.  1 

1.22 

6 

<2.l',0) 

Aj,  A,.  E 

A„Aj.  E 

0,  1 

1.45 

38 

■'  This  refen  to  the  irreducible  representation  of  the  P,  permutation  group  to  which  the  electronuclear  wave  function,  excluding  the 
nuclear  spin  part,  belongs  (31,321. 

*•  With(oul)  GP  refers  to  the  case  in  which  the  efTect  of  the  geometrical  phase  associated  to  conical  intersection  between  the  two  lowest 
electronic  state  potential  energy  surface  is  (not)  included  (31.32). 
cl  When  two  values  of  /7 are  indicated,  the  first  (second)  one  is  associated  with  the  first  (second)  permutation  symmetry 


to  the  lower  accuracy  of  the  32-bit  arithmetic  of  the 
former  with  respect  to  the  64-bit  arithmetic  of  the 
latter. 

The  absolute  times  presented  in  table  2  and  3  are 
apt  to  decrease  as  the  codes  are  improved  and  the 
numerical  parameters  are  further  tuned.  As  a  result, 
they  are  not  well  suited  for  an  appropriate  compare 
ison  of  the  relative  effectiveness  of  different  reactive 
scattering  methodologies  [8-19],  The  relevant  in¬ 
formation  in  those  tables  is,  instead,  the  relative  times 
among  different  machines  as  given  by  the  corre¬ 
sponding  speeds.  These  are  indicative  of  the  relative 
effectiveness  of  these  machines  for  performing  the 
reactive  scattering  calculations  described  in  this 
paper. 

The  efficiency  («)  of  the  parallel  LHSF  code  was 
determined  using  the  definition  t^TJNTN,  where 
f,  and  Ts  are  respectively  the  implementation  times 
using  a  single  processor  and  iV  processors.  The  single 
processor  times  are  obtained  from  runs  performed 
after  removing  the  overhead  of  the  parallel  code,  i.e. 
after  removing  the  communication  calls  and  some 
logical  statements.  Perfect  efficiency  ( t  =  1 .0 )  im¬ 
plies  that  the  .V processor  hypercube  is  N  times  faster 
than  a  single  processor.  In  fig.  4  efficiencies  for  the 
surface  function  code  (including  the  calculation  of 


the  overlap  and  interaction  matrices)  as  a  function 
of  the  size  of  the  primitive  basis  set  are  plotted  for 
2,  4,  8,  16,  32  and  64  processor  configurations  of  the 
hypercube.  The  global  dimensions  of  the  matrices 
used  are  chosen  to  be  integer  multiples  of  the  num¬ 
ber  of  processor  rows  and  columns  in  order  to  insure 
load  balancing  among  the  processors.  Because  of  the 
limited  size  of  a  single  processor  memory,  the  effi¬ 
ciency  determination  is  limited  to  32  primitives.  As 
shown  in  fig.  4,  the  efficiencies  increase  monotoni- 
cally  and  approach  unity  asymptotically  as  the  size 
of  the  calculation  increases.  Converged  results  re¬ 
quire  large  enough  primitive  basis  sets  so  that  the  ef¬ 
ficiency  of  the  surface  function  code  is  estimated  to 
be  about  0.95  or  greater. 

The  data  for  the  logarithmic  derivative  code  given 
in  table  3  for  a  245  channel  (i.e.  LHSF)  example 
show  that  the  Mark  lllfp  has  a  speed  about  62%  to 
that  of  the  CRAY-2but  only  about  31%  of  that  of  the 
CRAY  X-MP/48.  This  code  is  dominated  by  matrix 
inversions,  which  are  done  in  optimized  assembly 
code  in  all  three  machines.  The  reason  for  the  slow¬ 
ness  of  the  hypercube  with  respect  to  the  CRaYs  is 
that  the  efficiency  of  the  parallel  logarithmic  deriv¬ 
ative  code  is  0.52.  This  relatively  low  value  is  due  to 
the  fact  that  matrix  inversions  require  a  significant 
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Table  2 

Performance  of  the  surface  function  code  *' 


J 

Mark  Illfp01 

64  processors 

CRAY  X-MP/48 

CRAY-2 

lime  (h) 

speed  ( M  flops ) 

lime  (h) 

speed  ( Mflops) 

lime  ( h ) 

speed  ( Mflops) 

0 

0.71  “ 

too  •*» 

0.74“ 

96° 

0.49“ 

145“ 

1 

2.88° 

112“ 

3.04'' 

106° 

2.01  “ 

160“ 

2 

5.60*' 

124“ 

5.94 -» 

117“ 

3.96“ 

176“ 

"  This  code  calculates  the  surface  functions  at  the  $  I  values  of  fi  from  2.0  to  1 2.0  bohr  in  steps  of  0.2  bohr,  the  corresponding  overlap 
matrices  between  consecutive  values  of  0  and  the  propagation  matrices  in  p  steps  of  0. 1  bohr.  The  number  of  primitives  used  for  each 
J  and  described  in  the  remaining  footnotes  permits  us  to  generate  enough  LHSF  to  achieve  the  accuracy  described  in  the  text. 

“  Sixty-four  single  precision  processors. 

"  For  80  A|,  80  A2  and  160  E  primitives.  This  basis  is  larger  than  the  one  described  in  e)  below  and  is  needed  to  generate  the  same 
number  of  linearly  independent  surface  functions  as  in  e).  The  reason  for  this  difference  is  the  32-bit  arithmetic  of  the  Mark  IHfp 
compared  to  the  64-bit  arithmetic  of  the  CRAY  X-MP/48. 

“  Estimated  on  the  basis  of  the  absolute  measured  speed  on  the  CRAY  X-MP/48  and  the  measured  relative  speeds  of  the  Mark  IHfp 
with  respect  to  the  CRAY  X-MP/48. 

“  For  76  A|,  76  A2  and  152  E  primitives. 

'  ’  Measured  using  the  hardware-performance  monitor  of  the  PERFMON  and  PERFPRT  subroutines. 

11  This  time,  for  the  same  primitives  as  describes  in  e )  was  estimated  on  the  basis  of  the  relative  speeds  of  the  CRAY-2  and  CRAY  X- 
MP/48  measured  for  a  set  of  five  values  of  0.  It  is  smaller  than  the  time  in  e)  for  the  reason  in  h ). 

Estimated  on  the  basis  of  the  relative  speed  of  the  CRAY-2  with  respect  to  the  CRAY  X-MP/48  described  in  g).  The  reason  this  speed 
is  j  of  the  corresponding  CRAY  X-MP/48  speed  is  that  the  dominant  parts  of  the  calculation  are  optimized  assembly  code  matrix- 
vector  multiplications  for  which  the  CRAY-2  is  50%  faster  than  the  CRAY  X-MP/48.  Otherwise,  the  CRAY-2  is  slightly  slower  than 
CRAY  X-MP/48.  See  text. 

”  For  72  A|.  80  A2  and  152  E  primitives  of  even  parity  and  152  A,,  160  A2  and  312  E  primitives  of  odd  parity  These  numbers  of 
primitives  are  larger  than  the  ones  given  in  j )  for  the  reason  given  in  c ). 

11  For  64  a,,  76  A2  and  140  E  primitives  of  even  parity  and  140  Au  152  A2  and  292  E  primitives  of  odd  parity. 

“  Estimated  on  the  basis  of  the  relative  speeds  of  the  CRAY  X-MP/48  and  CRAY-2  and  the  measured  CRAY  X-MP/48  times  or  speeds. 

"  For  216  A,.  232  A2  and  448  E  primitives  of  even  parity  and  136  A,,  152  A2  and  288  E  primitives  of  odd  parity.  These  numbers  are 
larger  than  those  in  o )  for  the  reason  given  in  c). 

“'This  time  is  estimated  as  in  k).  since  the  calculation  cannot  be  done  on  the  CRAY  X-MP/48  because  of  insufficient  memory. 

*'  Estimated  to  be  the  same  as  in  f)  since  the  calculation  cannot  be  done  on  the  CRAY  X-MP/48  for  the  reason  given  in  m ). 

*'  For  204  A„  2 16  A2  and  420  E  primitives  of  even  parity  and  1 28  At,  140  A2  and  268  E  primitives  of  odd  parity. 


amount  of  inter-processor  communication.  Fig.  5 
displays  efficiencies  of  the  logarithmic  derivative 
code  as  a  function  of  the  number  of  channels  prop¬ 
agated  for  different  processor  configurations,  as  done 
previously  for  the  Mark  III  [25,42  J  hypercubes.  The 
data  can  be  fit  well  be  an  operations  count  formula 
developed  previously  for  the  matrix  inversion  pan 
of  the  code  [43];  this  formula  can  be  used  to  ex¬ 
trapolate  the  data  to  larger  numbers  of  processors  or 
larger  numbers  of  channels.  It  can  be  seen  that  for  an 
8  processor  configuration,  the  code  runs  with  an  ef¬ 
ficiency  of  0.81.  This  observation  suggested  that  we 
divide  the  Mark  IHfp  into  8  clusters  of  8  processors 
each  and  perform  calculations  for  different  energies 
in  different  clusters.  The  corresponding  timing  in¬ 
formation  is  also  given  in  tabic  3.  As  can  be  seen  from 
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the  last  row  of  this  table,  the  speed  of  the  logarithmic 
derivative  code  using  this  configuration  of  the  64 
processor  Mark  IHfp  is  48.5  Mflops,  which  is  about 
44%  of  that  of  the  CRAY  X-MP/48  and  88%  of  that 
of  the  CRAY-2.  As  the  number  of  channels  in¬ 
creases,  the  numbei1  of  processors  per  cluster  may  be 
made  larger  in  order  to  increase  the  amount  of  mem¬ 
ory  available  in  each  cluster.  The  corresponding  ef¬ 
ficiency  should  continue  to  be  adequate  due  to  the 
larger  matrix  dimensions  involved. 

In  the  near  future,  the  number  of  processors  of  the 
Mark  IHfp  will  be  increased  to  128  and  the  I/O  sys¬ 
tem  will  be  replaced  by  high  performance  CIO  (con¬ 
current  I/O)  hardware.  The  new  Weitek  coproces¬ 
sors  installed  since  the  present  calculations  were  done 
perform  64  bit  floating  point  arithmetic  at  about  the 
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Table  3 

Performance  of  tne  logarithmic  derivative  code  *’ 


Mark  Illfp  *” 


CRAY  X-MP/48  CRAY-3 


64  processor  8  clusters  of 

global  configuration  “  8  processors  “  ’ 


total  time  (h) 

4.8  •’ 

3.4f'*’ 

1.5 

2.9  ■" 

time  for  1  energy  (min) 

2.2" 

1.6" 

0.7 

1.3 

efficiency 

0.52 

0.81 

- 

- 

speed’’  (Mfiops) 

34.4  k’ 

48.5 k) 

110 

55.4 

*’  Based  on  a  calculation  using  245  surface  functions  and  131  energies,  and  a  logarithmic  derivative  integration  step  of  0.01  bohr. 

*'  Sixty-four  single  precision  processors. 

c>  The  calculation  for  each  energy  was  distributed  among  all  64  processors. 

'  ’  The  hypercube  was  configured  into  8  clusters  of  8  processors  each.  Each  cluster  did  full  calculations  for  1 6  energies,  for  a  total  of  1 28 
energies.  The  limes  reported  were  multiplied  by  1 3 1  / 1 28  for  normalization  purposes.  All  8  clusters  operated  simultaneously. 

*’  This  includes  1 .9  h  of  I/O  time. 

°  This  includes  1.6  h  of  I/O  time.  This  time  is  shorter  than  that  in  e)  because  of  a  different  and  more  efficient  broadcast  of  the  data 
between  the  host  and  the  8  clusters. 

"  Each  cluster  did  full  calculations  for  16  energies  for  a  total  of  128  energies.  The  total  time  reported  was  obtained  by  subtracting  the 
I/O  time  from  the  measured  time,  multipling  the  result  by  131/1 28  for  normalization  to  1 31  energies  and  adding  the  I/O  time. 

1,1  Estimated  on  the  basis  of  the  CRAY  X-MP/48  times  and  the  ratio  of  the  speeds  of  the  CRAY-2  and  CR  AY  X-MP/48  for  the  logarith¬ 
mic  derivative  code. 

"  This  includes  the  pro-rated  I/O  contribution. 

11  All  speeds  include  I/O  contribution. 

kl  Estimated  on  the  basis  of  the  measured  CRAY  X-MP/48  speed  for  the  logarithmic  derivative  code  and  the  relative  speeds  of  the  Mark 
Illfp  and  CRAY  X-MP/48  for  this  code. 


Global  Matrix  Dimension 


Fig  4.  Efficiency  of  the  surface  function  code  (including  the  calculation  of  the  overlap  and  interaction  matrices)  as  a  function  of  the 
global  matrix  dimension  ( i.e.  the  size  of  the  primitive  basis  set )  for  2,  4.  8,  16,  32,  and  64  processors.  The  solid  curves  are  straight  line 
segments  connecting  the  data  points  for  a  fixed  number  of  processors  and  are  provided  as  an  aid  to  examine  the  trends. 
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Fig.  5.  Efficiency  of  logarithmic  derivative  code  as  a  function  of  the  global  matrix  dimension  ( i.e.  the  number  of  channels  of  LHSF )  for 
8.  16,  32.  and  64  processors.  The  solid  curves  are  straight  line  segments  connecting  the  data  points  for  a  fixed  number  of  processors  and 
are  provided  as  an  aid  to  examine  the  trends. 


Table  4 

Overall  speed  of  reactive  scattering  codes  on  several  machines 


Mark  lllfp 

64  processor 

1 28  processors 

CRAY  X-MP/48 

CRAY -2 

CRAY  Y-MP/864 

surface  function  code 
for  7*2  (Mfiops) 

124 

240" 

117" 

176" 

232" 

logarithmic  derivative 
code"  (Mfiops) 

48.5  41 

127 

110s' 

55.4" 

187" 

total  main  memory  of 
computer  (64-bit  M words) 

32 

64 

8 

256 

64 

*'  Estimated  on  the  basis  of  the  64  processor  performance.  • 

"  For  single  processor  operation. 

"  For  243  channels  As  the  number  of  channels  increases,  the  Mark  lllfp  speed  increases  by  a  factor  not  excedmg  1.25.  but  the  speed  of 
the  CRAY  machines  remains  approximately  constant. 

1,1  Hypercube  configured  in  clusters  of  8  processors. 

"  This  speed  assumes  four-fold  increase  in  the  I/O  data  rate,  compared  to  the  64  processor  machine,  due  to  concurrent  I/O  hardware. 


same  nominal  peak  speed  as  the  32  bit  boards.  From 
the  data  in  the  present  paper  it  is  possible  to  predict 
with  good  reliability  the  performance  of  this  up¬ 
graded  version  of  the  Mark  lllfp.  A  CRAY  Y-MP/ 
864  has  just  been  installed  at  the  San  Diego  Super¬ 


computer  Center.  Initial  speed  measurements  show 
that  it  is  2  times  faster  than  the  CRAY  X-MP/48  for 
the  surface  function  code  and  1.7  times  faster  for  the 
logarithmic  derivative  code.  In  table  4,  we  summa¬ 
rize  the  available  or  predicted  speed  information  for 
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the  present  codes  for  the  current  64  processor  and 
near  future  128  processor  Mark  Illfp  as  well  as  the 
CRAY  X-MP/48.  CRAY-2  and  CRAY  Y-MP/864 
supercomputers.  Its  can  be  seen  that  Mark  Illfp 
machines  are  competitive  with  all  of  the  currently 
available  CRAYs  (operating  as  single  processor 
machines). 


5.  Summary 

We  performed  quantum  mechanical  reactive  scat¬ 
tering  calculations  on  the  Mark  Illfp  hypercube  par¬ 
allel  computer.  The  results  obtained  for  the  H  +  H2 
system  /=0,  1,2  partial  waves  agree  well  with  those 
from  a  CRAY  X-MP/48  and  a  CRAY-2.  The  reso¬ 
nance  structure  in  the  J=1  calculations  is  consistent 
with  a  selection  rule  developed  previously  (9,16]. 
The  high  degree  of  parallelism  of  the  most  time-con¬ 
suming  step  of  the  surface  function  calculation  (the 
evaluation  of  two-dimensional  numerical  quadra¬ 
tures)  leads  to  a  high  efficiency  for  that  calculation. 
As  a  result,  the  speed  of  the  64  processor  Mark  Illfp 
for  the  surface  function  calculation  is  about  the  same 
as  that  of  the  CRAY  X-MP/48  and  about  0.7  of  that 
of  the  CRAY-2.  When  configuring  the  Mark  Illfp  into 
8  clusters  of  8  processors  each,  the  logarithmic  de¬ 
rivative  code  is  about  56%  slower  than  the  CRAY  X- 
MP/48  and  12%  slower  than  the  CRAY-2.  The  speed 
of  the  128  processor  Mark  Illfp  soon  to  become 
available  should  exceed,  both  for  the  surface  func¬ 
tion  calculation  and  the  logarithmic  derivative  cal¬ 
culation,  those  of  the  CRAY  X-MP/48  and  CRAY- 
2;  however,  although  still  comparable  to  the  CRAY 
Y-MP/864  for  the  surface  function  code,  it  will  be 
32%  slower  for  the  logarithmic  derivative  code  (the 
CRAYs  operating  as  single  processor  machines). 
These  results  demonstrate  the  feasibility  of  perform¬ 
ing  reactive  scattering  calculations  with  high  effi¬ 
ciency  in  parallel  fashion.  As  the  number  of  proces¬ 
sors  continues  to  increase,  such  parallel  calculations 
in  systems  of  greater  complexity  will  become  prac¬ 
tical  in  the  not  too  distant  future. 
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